The Inclined Wooding Problem

A constant-thickness thermal boundary layer is formed by the presence of uniform suction into a constant temperature hot surface bounding a porous medium—the Wooding problem. When the hot surface lies below the porous medium, the system is susceptible to thermoconvective instability. The present paper is concerned with how the classical linear stability analysis is modified by inclining the heated surface. Analysis is confined to disturbances in the form of transverse rolls because the equivalent analysis for longitudinal rolls may be described analytically in terms of that for the horizontal layer. The linear stability analysis is made difficult by the absence of a second bounding surface, and the method of multiple shooting is needed in order to obviate the consequence of having a stiff system of disturbance equations. Therefore, the computational domain is split into 5 or 10 subdomains. It is found that all modes of instability travel up the heated surface unless the surface is horizontal. The system is found to be linearly stable for all inclinations above 31.85473∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$31.85473^\circ $$\end{document}, a value which is remarkably close to that for the inclined Darcy-Bénard problem.


Introduction
The Wooding problem [which was named as such by Pieters and Schuttelaars (2008)] may be described as one of the classical thermoconvective instability problems in fluid mechanics, one which is of Rayleigh-Bénard type, a thermoconvective instability. Wooding noted how the upward ground water movement which had been observed in Wairakei, New Zealand, when combined with surface cooling and evaporation, produces conditions whereby relatively dense fluid lies above relatively light fluid, and hence this provides the potential for thermoconvective instability to occur (Wooding 1960). The idealised problem shares some features with the Darcy-Bénard problem, but the two chief differences are the background vertical fluid motion and the fact that the porous region is of semi-infinite extent, rather than being of finite extent. The background fluid motion causes a thermal boundary layer to form naturally, and it is this which is potentially unstable. The semi-infinite nature of the domain means that convective motions may be found well outside of the thermal boundary layer. Sutton (1970) considered the effect of a weak throughflow on the Darcy-Bénard problem in a finite box, and she provided numerical and series solutions for the variation of the critical Darcy-Rayleigh number with the aspect ratio of the box. In a significant paper, Homsy and Sherwood (1976) extended Sutton's work numerically to a wide range of throughflow Péclet numbers and also provided the first energy stability analysis for this type of flow. They found that the linear and energy analyses gave the same critical value of Ra, the Darcy-Rayleigh number, only when there is no throughflow; otherwise, the value of Ra c for the energy analysis is below that for a linear analysis, implying that subcritical convection arises due to the presence of throughflow. Pieters and Schuttelaars (2008) presented a very detailed account of nonlinear flow for Sutton's configuration for a range of Péclet numbers. Using a projection method, detailed bifurcation diagrams are created to determine the possible range of flows. They also derive a weakly nonlinear theory, and their Fig. 2 suggests that the transition from supercritical onset to subcritical onset takes place when the Péclet number is roughly 4.1. An alternative weakly nonlinear theory may be found in Riahi (1989) where the heated surface is assumed to have an almost perfectly constant heat flux boundary condition. This mathematical device allowed the author to develop a three-dimensional analysis, and he was able to determine that both up-hexagons and square planforms are possible convection patterns.
Many other papers have appeared on this general topic, but these have provided information on how other effects affect the onset criteria, rather than being devoted to furthering the essential pure Wooding and Sutton problems beyond what has been presented above. Such papers include the following. Different boundary conditions were considered by Nield (1986). Shivakumara (1999) added the Brinkman and the advective inertia terms to Darcy's law, while Khalili and Shivakumara (2003) also included the Forchheimer terms. Khalili and Shivakumara (1998), Yoon et al. (1998) and Deepika et al. (2017) considered internal heat generation within the layer, while Nield andKuznetsov (2013, 2015) allowed layering. Shivakumara and Khalili (2001) considered double diffusion, i.e. two diffusing components, while Shivakumara andSureshkumar (2007, 2008), Barletta and Storesletten (2016) and Rees (2018) analysed some non-Newtonian fluids. Hill (2007) has used an energy stability analysis to determine the effect of cubic Forchheimer terms on convection, while Hill et al. (2007) considered the effect of the water density maximum at 4 • C. Rees (2009) considered a forced convection boundary layer flow in the presence of surface suction, and even in this case it is possible to obtain subcritical vortex convection given a sufficiently large upstream disturbance. Local thermal nonequilibrium effects were found by Patil and Rees (2013) to cause double minima in the neutral curve under some circumstances. Anisotropic effects were considered by Rees and Storesletten (2002). Solute immobilisation effects, such as one would find in mineral trapping for CO 2 sequestration, have been considered by Maryshev and Lyubimova (2016). Ferroconvection was studied by Nanjundappa et al. (2015).
The present paper is concerned with the effects of inclination on the Wooding problem, and attention is confined to transverse rolls. This type of variation is not without its precedents, for Rees and Bassom (2000), Rees and Postelnicu (2001), Rees and Storesletten (2002) Rees and  and Barletta and Storesletten (2011) all consider variants on the inclined Darcy-Bénard problem. Unlike the corresponding clear-fluid analogues, convection in the form of transverse rolls appears only within a certain range of inclinations about the horizontal. Indeed, Gill (1969), Lewis et al. (1995) and Straughan (1998), using different methods, have all shown that all disturbances decay when the layer is vertical. More specifically, it appears to be a general property of inclined porous layers that convection is extinguished once a critical inclination is reached [although this may not be true for hydrostatic pressure boundary conditions: see Barletta (2015)]. The overall aim for the present paper is to compare the stability properties of a semi-infinite system with those for finite-depth layers. Specifically we shall present the neutral curves, determine the variation of the critical values with inclination and find the isola point at which all convection ceases.

Governing Equations and Basic State
We are considering the stability of the convective flow which is induced by the presence of an inclined heated plate with uniform surface suction and which is embedded within a homogeneous isotropic porous medium. The full three-dimensional governing equations are given by the equation of continuity, Darcy's law as modified by buoyancy forces subject to the Boussinesq approximation, and the equation of heat transport. When in dimensional form these respective equations are, Although all the terms introduced above take their common meanings for this type of context and are listed in the Nomenclature, we note specifically that the heated surface is aligned at an angle, α, to the horizontal, that the x and z axes are aligned parallel with and perpendicular to the heated surface, respectively, and that the y axis is horizontal. The value, σ , is the ratio of the heat capacity of the saturated medium to that of the fluid: where φ is the porosity. The boundary conditions are that, In the above, we need to have W > 0 so that the heated surface is a suction surface.
The nondimensionalisation of Eqs.
(1) to (5) is achieved using the following substitutions, where is a lengthscale based upon the thermal diffusivity and the suction velocity; this value is also equivalent to the width of the resulting thermal boundary layer, and therefore it forms a natural lengthscale. Hence, we obtain the non-dimensional equations, The quantity, is the Darcy-Rayleigh number, where L is given by (10). The boundary conditions now become, The basic state which we will analyse for stability is obtained by first assuming that the flow is steady and dependent only upon z. Therefore, we obtain, This paper will be devoted to a two-dimensional linear stability analysis where disturbances take the form of transverse rolls, and this may be analysed by the introduction of a suitable streamfunction. Therefore, we set, and so Eqs. (11) to (13) reduce to the forms, Within this formulation, the basic state is given by, We note that the constant term in the expression for ψ has been chosen to be such that ψ = 0 on z = 0, but it could take any constant value. We note also that the basic temperature field 123 has an O(1) thickness which is independent of the angle of inclination of the heated surface.
On the other hand, the magnitude of the induced flow field, which is parallel to the heated surface, depends strongly on α. Its strength is maximised when the surface is vertical, and no flow is induced (apart from the suction velocity) when it is horizontal.

Linear Stability Analysis
We now proceed to describe the linear stability of the basic state to two-dimensional disturbances. To that end we shall substitute, into Eqs. (18) and (19). The disturbance quantities, and , are assumed to have asymptotically small amplitudes, and therefore products of the disturbances are neglected. The resulting linearised stability equations are now, where the boundary conditions are that both and are zero on the heated surface and tend to zero far from the surface. These equations may now be split into independent monochromatic modes in the x-direction by setting, where k is a wavenumber and λ = λ r + iλ i is the complex exponential growth rate for disturbances. The complex equations for f and g take the forms, and the boundary conditions are

General Observations
Equations (25)-(27) may be regarded either as (i) an ordinary differential eigenvalue problem for the complex growth rate, λ, as a function of Ra, k and α, or else (ii) as an ordinary differential eigenvalue problem for the Rayleigh number and λ i as a function of k and α by setting λ r = 0 to marginally stabile conditions. The former point of view lends itself naturally to the solution procedure that was adopted by Rees and Bassom (2000) whereby all derivatives are approximated using second-order accurate finite differences, and the resulting ordinary differential eigenvalue problem is thereby reduced to that of a matrix eigenvalue problem. Then values of λ are computed on a fine grid of values of Ra and k, and the neutral curves are obtained by contouring where λ r = 0. This has the advantage that all neutral curves, including any isolated loops, are obtained. When taking the latter point of view, the equations may be solved by more conventional methods, such as a Runge-Kutta scheme and a shooting methodology in order to find Ra in terms of the other parameters. It was found that the general properties of these equations made it substantially more difficult to obtain accurate solutions than for the inclined Darcy-Bénard layer which was considered in Rees and Bassom (2000). There are a variety of reasons for this. Concentrating first on the matrix eigenvalue methodology, the fact that the domain is semi-infinite (as opposed to being of uniform thickness, as in Rees and Bassom 2000) means that a substantially larger number of grid points needs to be employed, and this leads to a very much slower computation. In addition, when k is large, then disturbances tend to occupy relatively narrow regions when compared with the basic state, and therefore the resolution of these regions becomes difficult. When k is small, then the decay rate for the streamfunction as z increases is rather small because the asymptotic profile for f may be shown to be proportional to e −kz . This combination of difficulties led us to take the more conventional Runge-Kutta/shootingmethod route.
In a sense, the same difficulties exist when using the fourth-order Runge-Kutta method, except that solutions were obtained very much more quickly. The slowly decaying behaviour of f when k is small was modelled by replacing the f → 0 as z → ∞ boundary condition with the Robin condition, The presence of the g term in Eq. (26) means that the e-folding distance for the temperature disturbances is of the order of unity, and therefore the value, z max = 10, for example, means that g has often decayed to below 10 −4 in magnitude. Therefore, the boundary condition in (28) accounts perfectly for the asymptotic exponential decay behaviour of f even when k is extremely small. When k takes large values, the resolution of the narrow near-surface regions is not merely a matter of taking a sufficiently large number of intervals. However, the numerical problem now becomes stiff, and therefore our chosen method for overcoming this last difficulty was to adopt the variant known as multiple shooting. In the standard shooting method, Eqs. (25) and (26) form a 4th-order complex system, or equivalently, an 8th-order real system. The addition of two extra equations, Ra = 0 and λ i = 0, renders us with a 10th-order system overall. The adoption of the multiple-shooting technique was generally undertaken with five subregions being solved simultaneously, and with appropriate interface conditions being imposed-details are given below. This yields a system of 42nd-order with 37 unknown initial conditions, but despite its size it was found to yield accurate solutions quite quickly. Of course, even this complexity only delays to larger values of k the recurrence of the problem of stiffness, but we found that 5 subregions were sufficient for almost all our needs, and some computations with 10 subregions were undertaken to check that the methodology had been applied correctly.

Numerical Methodology
Our chosen basic method is the shooting method where the solver is the classical fourthorder Runge-Kutta method. Given the above discussion, no consistent value of z max was used, but this depended to some extent on which part of parameter space (Ra, k and α) was being occupied, but at the earliest stages tests were conducted to assess how large the 123 computational domain should be and how many grid points were to be employed. Generally z max was either 8 or 10, and further increases produced insignificant changes in the solution. With regard to the steplengths, our solutions are correct to five significant figures or better.
Neutral curves, which correspond to nonzero solutions to Eqs. (25) and (26), are found by setting λ r = 0. These are found not always to be single-valued functions of the wavenumber, but can exhibit vertical tangents and therefore we have adopted a curve-following (or pseudoarc-length) strategy to trace out the curves. At its most straightforward this means that we need to solve an 11th-order system of equations, as follows.
where Ra and k will be replaced by v 9 and v 11 , respectively, in the numerical codes but are retained here for clarity. The above equations were solved subject to the boundary conditions, ⎛ The entries for z = 0 which are marked with a black square are those which are unknown and which need to be found by insisting that the given conditions at z = z max are satisfied. The values for v 4 and v 8 at z = 0 (i.e. g (0) = 1) are the normalisation conditions to ensure nonzero solutions, while those involving k at z = z max are equivalent to Eqs. (28). Finally, the quantity, S, is defined as, In this expression k old and Ra old are the values of k and Ra at the previously computed point on the neutral curve, and the values δ k and δ R are the largest increments in k and Ra allowed when moving along the neutral curve; typical values are δ k = 0.001 and δ R = 0.5. Satisfaction of this condition allows the code to compute the neutral curve past turning points. The most important point on a neutral stability curve is its minimum. These points may be found by solving the system, Eqs. (25) and (26), together with that system which is formed by taking its k-derivative where ∂Ra/∂k = 0 is set. So if we set, F = ∂ f /∂k and G = ∂ g/∂k, then this extra system is given by, and the boundary conditions are In the above = ∂λ/dk. When minima are sought, Eqs. (32) and (33) are also converted into first-order form and solved simultaneously with Eq. (29). When using the standard shooting method, this extended system forms one of 20th-order with eight unknown initial conditions. Solution of this system yields Ra c , k c and the values of λ and ∂λ/∂k as functions of α. The present stability problem admits both an isola point and a saddle point and these points require that both ∂Ra/∂k = 0 and ∂Ra/∂α = 0 are satisfied simultaneously at isolated points in parameter space. The appropriate system which needs to be added to both Eqs. (25) and (26) and Eqs. (32) and (33) may be obtained by differentiating Eqs. (25) and (26) with respect to α and by setting ∂Ra/∂α = 0. If we set F = ∂ f /∂α and G = ∂ g/∂α, then we have, subject to the boundary conditions, where L = ∂λ/∂α. Finally, for large values of k and for large values of Ra the equations become too stiff to be able to be solved in a sufficiently large computational domain that good accuracy is preserved, and therefore we resorted to using the method of multiple shooting using five subdomains, as mentioned above. As an example of this, we could define five equal domains lying between z = 0 and z = 10 and name them subdomains 1 through to 5, where 0 ≤ z ≤ 2 in subdomain 1. The first eight variables shown in Eq. (29) have their counterparts in all the subdomains. For ease of bookkeeping, subdomain 1 contains variables 1 to 8, subdomain 2 contains variables 9 to 16 and so on. These 40 variables, together with the three corresponding 123 to Ra, k and λ i , which are designated variables 41 to 43, are then solved in the region 0 ≤ z ≤ 2 simultaneously, where appropriate adjustments are made to the exponential functions, e −z , in domains 2 to 5. The initial conditions given in Eq. (30) now apply in subdomain 1 and the final conditions given there apply in subdomain 5. Then it is necessary to apply continuity conditions between the subdomains; if we were to consider the interface between subdomains 1 and 2, then the continuity conditions would be, In this case, we would have 43 equations with 38 unknown initial conditions. A similar multiple-shooting idea was used for determining the minima (and, in some cases, maxima) in the neutral curves, but without the curve-tracking idea being implemented. In this case, we have 84 equations with 72 unknown initial conditions, and when we adopted 10 subdomains, we had 164 equations with 152 unknown initial conditions.

Neutral Curves
A selection of neutral curves is presented in Fig. 1 for a chosen set of values of the inclination, α. When α = 0 we reproduce the Wooding problem which has been studied extensively by authors such as van Duijn et al. (2002) and Pieters and Schuttelaars (2008). When the heated surface is horizontal (α = 0) The minimum of the neutral curve corresponds to Ra c = 14.35219, k c = 0.75887; which upgrades slightly the data given in Rees (2009). As the inclination increases, the critical Rayleigh number increases because the vertical component of the buoyancy force diminishes. The critical wavenumber also increases and does so almost linearly at first. The narrow black line in Fig. 1 which marks the largest inclination at which one has incipient convective instability; it is the lower of the two black disks shown in Fig. 1. Such an isola point is a general feature of inclined thermoconvective stability problems in porous media, and Table 1 reproduces details of these points for different types of inclined layers.
Here CT refers to constant temperature boundary conditions, CHF to constant heat flux and MIX to one of one type and one of the other. In addition, DB refers to a Darcy-Bénard-like heating from below, while IHG refers to internal heat generation within the layer. The very obvious closeness of the maximum inclinations shown in Table 1 appears to suggest that there might be a deep connection between these five convective instability problems. However, the study of an anisotropic version of the CT-DB problem by Rees and Postelnicu (2001) shows that α max can achieve values as low as 12 • and as high as 47 • depending on the degree of anisotropy of the porous medium. Table 1 also shows that there are some significant variations in the values of Ra c and k c between the different problems represented. All three of those labelled, DB, are layers with unit thickness and heated from below, and therefore there isn't a large variation in the Red curves correspond to α = 21 • to 31 • in steps of 1 • . Continuous blue curves correspond to α = 29.2 • to 29.8 • in steps of 0.2 • , while the dashed blue curve is α = 29.3 • . The dashed green curve corresponds to 31.5 • . The black circles denote the isola and saddle points, while the black lines passing through them are the loci of extrema in the neural curves which are also labelled according to whether they form the main, left or right-hand branches  Barletta et al. (2014) respective values of Ra c and k c . The two labelled, IHG, correspond to situations where the temperature gradient is destabilising over only a lower fraction of the layer; thus Ra c must be larger and, given that cells will occupy mainly those regions where there is a negative temperature gradient, the wavenumbers will be substantially larger. For the present Wooding problem, the basic temperature field occupies a region which is larger than unity, and therefore both Ra c and k c will be smaller than for the DB problems. While a superficial glance at the general shapes of the neutral curves given here and those shown in Rees and Bassom (2000) might suggest that they are very similar, there are nevertheless some very important differences. One of these is that every neutral curve shown here (apart from when the heated surface is horizontal) corresponds to an onset mode which travels up the heated surface. This is due to the lack of symmetry in the basic state, a property which is shared with the configurations studied in Barletta and Storesletten (2011). The computed values of λ i are always negative here, and therefore the speed of any convection pattern in the x-direction, which is given by −λ i /k [see Eq. (24)], is always positive. On the other hand, the inclined Darcy-Bénard problem continues to exhibit stationary convection at each neutral curve minimum. Pairs of travelling modes do exist, but these bifurcate from the stationary branches at larger values of Ra, and therefore do not form the most unstable mode; a similar statement may be made for the CHF-DB problem studied by Rees and Barletta (2011).
The second difference is that the locus of extrema here consists of two disconnected branches, as opposed to the single branch shown in Rees and Bassom (2000). However, both configurations have a saddle point, the present one being located at and this is shown as the upper black disk in Fig. 1.

Variations with Inclination
Further detail which may be gleaned from Fig. 1 Figure 2 shows three branches and these correspond to extrema in the neutral curves in Fig. 1. The branches are labelled Main, Left and Right for ease of references. The Main branch shown in Fig. 2 is even about α = 0 and therefore negative inclinations are not presented. For each value of α on the main branch, linear instability corresponds to the region either above the main branch or, when α is sufficiently large, to that region between the lower and upper sections of the main branch. The right and left branches form part of the solution structure for the linearised eigenvalue probem, but as they lie above the lower part of the main branch, they play no role for the physical problem in reality, for they correspond to large values of Ra. The turning point in the main branch corresponds to the isola point.
The variation of k c with α, as shown in Fig. 3, is very modest when α increases from zero to the isola point (i.e. the turning point of the main branch shown here). To two decimal places this variation lies between 0.76 when the surface is horizontal and 1.14 when it is at an inclination of almost 30 • . In turn this means that the wavelength of the most unstable mode decreases by about 2 /3 from roughly 8.27 to 5.51; a wavelength decrease is a common qualitative feature of the effect of inclination on the onset of convection in porous media.   Figure 4 displays the variation with inclination of the relative speed of the cells, U, given in Eq. (42). The velocity of the cells up the layer has already been stated to be precisely, −λ i /k, but Eq. (16) shows that the maximum value of the basic velocity profile, which occurs when z = 0, takes the value, Ra sin α. Therefore, U represents the velocity of the cells compared with the maximum for the basic state. It will generally be the case that the larger is U, the more the disturbance is confined to the near-surface region.
Although the speed of the disturbances along the layer is zero when the heated surface is horizontal, the numerical data suggest that there is a well-defined limit, as α → 0, which is 0.2417. Thus, at very slight inclinations, the cells travel along the layer at roughly a quarter of the maximum velocity of the basic flow. As the inclination increases, both the wavenumber (see Fig. 3) and U increase. Thus U increases from 0.2417 to roughly 0.469 at the isola point. It is well-known that, for convective instability problems of boundary layer flows, the disturbance tends to be increasingly confined to the neighbourhood of the bounding surface when the wavenumber is large, and this may be traced mathematically to the e-folding distance of equations such as Eq. (25). Therefore, the disturbance tends to occupy regions where the background flow is faster, and hence the speed of the cells increases.

Disturbance Profiles
Finally, we display some disturbance profiles in Fig. 5. Four cases have been chosen, and these correspond to the respective minima in the neutral curves for α = 0 • , 20 • , 25 • and 31.5 • , the last one being very close to the isola point. Disturbance streamlines and isotherms have been created by taking the real parts of the functions given in Eq. (24). In each of the four diagrams composing Fig. 5 the height corresponds to z max = 10, while the horizontal width comprises precisely one period of the critical disturbance.  As one would expect, the streamfunction and the isotherms are exactly out of phase when the inclination is zero, for then there is no basic flow in the horizontal direction. As the heated surface is inclined increasingly, the strength of the basic flow also increases and this, together with the z-variation of that basic flow, has the effect of deforming the convection cells. We also see that the isotherms are confined to a narrower region close to the heated surface than are the streamlines, as discussed earlier. Generally, it is the case that the isotherms are confined to a region which is, at most, given by what is displayed for α = 0 • . When the wavenumber is large, it is possible for that region to be more restricted (such as for α = 31.5 • here) but it cannot be larger. The mathematical reason for the upper restriction on the width of the thermal disturbances is the presence of the g term in Eq. (26), whose presence is a direct consequence of surface suction. On the other hand, for those cases where k is quite small (not depicted here for reasons of brevity), the streamfunction field expands substantially, and this is why we adopted the boundary conditions given in Eq. (28) which accounts for the correct exponential decay in the isothermal region well outside where thermal effects are significant.

Conclusions
In this paper, we have sought to provide a comprehensive account of the effect of inclination on the classical Wooding problem. While other authors have also considered the onset of convection in similarly inclined configurations, the chief difference here is that there is no second bounding surface, and thus the flow and instabilities take place in what is effectively a deep-pool situation. Practically, this means that very considerable care had to be taken to maintain accuracy, to avoid difficulties due to the governing equations being too stiff to use the usual shooting method, and to account for slowly decaying streamfunction disturbances. The combined modifications to our numerical method meant that overall numerical scheme was robust.
The inclined Wooding problem displays some characteristics which are common to very many other inclined convective instability problems, including a maximum inclination above which the system is linearly stable. But there are other characteristics which are not the same, such as whether the disturbances are stationary modes or travelling modes. It was also seen that the symmetry of the basic state allows one to predict which property is assumed by the system.
The authors are well aware that the Wooding problem is a classic example of one where the bifurcation to strongly convecting flow is subcritical at the critical wavenumber (see Rees 2009). By continuity, this will also be true for at least a range of inclinations near to the horizontal. It is intended that the present analysis will be extended in order to determine if this remains true for all inclinations for which instability is possible, or whether there is an inclination above which the onset of convection is supercritical. This could be undertaken using a weakly nonlinear theory, energy methods or by solving the full nonlinear equations directly. It is hoped to report on this aspect in due course.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.