Horizontal Cellular Oscillations Caused by Time-Periodic Resonant Thermal Forcing in Weakly Nonlinear Darcy-Bénard Convection

: The onset of Rayleigh-Bénard convection in a horizontally unbounded saturated porous medium is considered. Particular attention is given to the stability of weakly nonlinear convection between two plane horizontal surfaces heated from below. The primary aim is to study the effects on postcritical convection of having small amplitude time-periodic resonant thermal forcing. Amplitude equations are derived using a weakly nonlinear theory and they are solved in order to understand how the ﬂow evolves with changes in the Darcy-Rayleigh number and the forcing frequency. When convection is stationary in space, it is found to consist of one of two different types depending on its location in parameter space: either a convection pattern where each cell rotates in the same way for all time with a periodic variation in amplitude (Type I) or a pattern where each cell changes direction twice within each forcing period (Type II). Asymptotic analyses are also performed (i) to understand the transition between convection of types I and II; (ii) for large oscillation frequencies and (iii) for small oscillation frequencies. In a large part of parameter space the preferred pattern of convection when the layer is unbounded horizontally is then shown to be one where the cells oscillate horizontally—this is a novel form of pattern selection for Darcy-Bénard convection.


Introduction
The classical Bénard problem as applied to a porous medium was first studied by Horton and Rogers (1945) [1] and independently by Lapwood (1948) [2]; these authors considered the onset of convection in a horizontal saturated porous layer heated uniformly from below.Using a linear stability theory, the neutral curve which relates the Darcy-Rayleigh number, Ra, to the wavenumber, k, may be shown to take the form, (see Rees 2000, for example) [3] where the value of n corresponds to the number of rolls which are stacked above one another in the layer.This variation of Ra with k is shown in Figure 1 and it may be shown easily that the minimum/critical value is Ra c = 4π 2 , which occurs when k = π and n = 1.These values correspond to rolls with a square cross-section.Neutral curve for the classical Darcy-Bénard problem with n = 1; see Equation (1). the critical values are given by Ra c = 4π 2 and k c = π.Palm et al. (1972) [4] subsequently analysed the moderately supercritical flow using a series representation in terms of powers of (Ra − Ra c )/Ra.A very detailed numerical stability analysis was undertaken by Straus (1974) [5] who delineated the region in (Ra, k)-space within which steady rolls form the stable planform of convection.
The general topic of Darcy-Bénard convection continues to grow rapidly as further effects, realism and/or practical application are added to the original Horton-Rogers-Lapwood configuration [2,6].Reviews have been presented by Rees (2000) [3] and Tyvand (2002) [7], but the topic has continued to be the subject of substantial interest in the time since then, and therefore we would refer the reader to the latest edition of Nield and Bejan (2017) [8] for an up-to-date account of the topic.
In the present paper, we shall confine our interest to the role played by boundary imperfections, itself a well-researched subtopic.It is widely recognised that completely uniform and idealised boundary conditions are not easily achieveable in practice, and therefore authors such as Riahi (1983) [9], Saleh et al. (2011) [10], Motjabi and Rees (2011) [1] and Rees and Mojtabi (2011) [11] have concerned themselves with more realistic boundary conditions where the perfectly conducting bounding surfaces are replaced by conducting solids.The resulting stability properties are influenced profoundly by the nature of the solids used and their thickness.Thus it is possible for perfectly conducting boundary conditions on the outer surfaces to give rise to stability properties which are more closely associated with constant heat flux surfaces, and vice versa.It is also possible for roll solutions to be unstable in some circumstances, in which case convection with a square planform forms the preferred pattern (see Riahi 1993, Rees andMojtabi 2011) [11].
A little later, Mamou et al. (1996) [21] and Banu and Rees (2001) [22] applied thermal boundary variations of the form cos π(x − Ut), which correspond to travelling thermal waves at the critical spatial wavenumber, π.The former paper considered strongly supercritical convection using fully numerical methods, while the latter applied a weakly nonlinear theory in order to understand the subtle details which arise closer to the onset of convection.The effect of having a forcing wave speed, U, is to cause a competition between the tendency of the rolls to remain stationary, which is what happens naturally in a layer without such imperfections, and the tendency to follow the forcing wave.It was found that there is a catastrophic transition between the wave-following regime and the quasi-stationary regime in that the amplitude and horizontal velocity of the convection pattern change suddenly as the Darcy-Rayleigh number increases.This system was also found to display hysteresis, for the reverse transition which occurs as the Rayleigh number decreases, does so at a different value of the Rayleigh number.
The present paper is concerned with a physical configuration for which convection cells also have a tendency to move.For the Darcy-Bénard problem one of the easiest ways in which this may arise is when the fluid is also subject to a horizontal pressure gradient.Prats (1967) [23] showed that the convection cells then move with precisely the same velocity as the pressure-gradient-induced flow which exists in the absence of heating.Thus the dynamics of the ensuing convection are precisely the same as that for a layer without the background flow when it is considered in a frame of reference which is moving with the background flow.An interesting variation on this Darcy-Bénard-Prats problem was studied by Dufour and Néel (1996) [24] where a uniform horizontal flow is injected into a horizontally semi-infinite porous layer heated from below.In general the flow remains roughly uniform near to injection and there exists a spatial transition towards the time-periodic flow described by Prats (1967) [23].The length of the transition region depends on the magnitudes of the injection velocity, Q (a Péclet number) and the Darcy-Rayleigh number, increasing as Péclet number increases but decreasing as the Darcy-Rayleigh number increases; these phenomena are related to the phase diffusion speed of the convection cells.A much more recent paper by Rees and Mojtabi (2013) [1] shows that the presence of solid conducting boundaries serves to reduce the phase speed of the convecting cells.This reduction may be attributed to what might be called a thermal drag which arises due to the fact that the time-varying convection mode has to move through a stationary boundary of finite height.
Another mechanism which relies neither on an external pressure gradient nor on moving boundary conditions to induce phase drift is when convection takes place in a cylindrical container with a circular planform.Although no studies have yet been published on such flows in porous media, it is well-known that such a variant of the classical Bénard problem may sometimes exhibit a rotating spiral convection pattern even though all the boundary conditions are steady and uniform; see Plapp et al. (1998) [25].
In this paper we consider thermal forcing of the form cos ωt cos πx, and this may be regarded as the standing wave counterpart to the work of Banu and Rees (2001) [22].Although this thermal forcing is now stationary in space, we shall show later that it remains possible for rolls to exhibit a time-periodic horizontal motion.As shall be seen, the reason for this periodic motion is that the preferred phase of the convection cells depends (in a quasi-static sense) on the sign of the forcing amplitude, cos ωt, and therefore the phase of the cells is drawn towards two different fixed points during the two different parts of the forcing cycle.Our subsequent analyses involve a derivation of the weakly non-linear amplitude equations, together with both numerical simulations and asymptotic analyses of those amplitude equations.Attention is also focused on the transition between two different types of stationary motion, and on a twitching motion at high Rayleigh numbers where cells execute small-amplitude periodic movements in phase.

Governing Equations
Darcy's law gives the most basic equation governing fluid flows in a porous medium.It is the macroscopic law which relates the fluid flux to the applied pressure gradient.For buoyancy-affected flows of a Boussinesq fluid, where the porous medium is assumed to be homogeneous, nondeformable and isotropic, the dimensional equations which govern the fluid motion are, In these equations all the terms take their common meanings.The coordinates x and y are in the horizontal directions, while z is vertically upwards.T is the temperature of the saturated medium, T c is the mean upper (cold) boundary temperature, p is the pressure, K the permeability, µ the viscosity, g gravity, ρ f a reference density of the fluid, and β the coefficient of cubical expansion of the fluid.The quantities u, v and w, correspond to the fluid seepage velocities in the x, y and z directions, respectively.
The full equations are completed by the equation of continuity, ∂u ∂x and the heat transport equation, where κ is the thermal diffusivity of the saturated porous medium and, is the heat capacity ratio of the saturated medium to that of the fluid.The value, c, is the specific heat and ρ the density of the medium, with the subscripts f and s referring to the fluid and solid phases, respectively.The nondimensionalisation of Equations ( 2)-( 4) is achieved using the following substitutions, where d is the height of the layer.Here, T c and T h are the cold (upper) and hot (lower) mean temperatures of the horizontal bounding surfaces.Hence, the non-dimensional equations are, ∂u ∂x In the above, the parameter, is the Darcy-Rayleigh number; the usual clear-fluid Rayleigh number may be obtained by replacing the permeability, K, by d 2 .The Darcy-Rayleigh number expresses the balance between buoyancy and viscous forces.
Assuming that the flow is two dimensional, then v = 0 and all the derivatives with respect to y may be neglected.Hence, by introducing the stream function, ψ, using u = −ψ z and w = ψ x , not only is Equation (8) satisfied but Equations ( 9) and (10) become, Both of the equations above are to be solved analytically using a weakly nonlinear analysis.
The boundary surfaces are considered to be impermeable.The thermal boundary conditions comprise a combination of heating from below and a small-amplitude variation which is periodic in both time and space.Therefore the boundary conditions are where δ 1 is the amplitude of the thermal forcing and ω is its nondimensional temporal frequency.

Weakly Nonlinear Analysis
In this section we shall carry out a weakly nonlinear analysis in order to obtain a complex Landau equation for the amplitude of convection in the presence of the thermal boundary imperfections given in Equation ( 14).This will be undertaken using the scalings which were used in Rees and Riley (1989a) [13] and Banu and Rees (2002) [22], and which are suitable for those cases where the wavenumber of the imperfection takes the critical value.We follow the methodology introduced by Newell and Whitehead (1969) [26] for the Bénard problem.
Therefore we shall set the amplitude of convection to be of O( ) and the Darcy-Rayleigh number shall be within a distance of O( 2 ) of its critical value, 4π 2 .Given that the imperfection has the critical wavenumber, we need to set δ = 3  1; see Rees andRiley (1986,1989a) [12,13].The appropriate time scale is of O( −2 ).Equations ( 12) and ( 13) may now be expanded using the following asymptotic series, and, given the above observations, we set where R 0 = 4π 2 .The leading term in Equation ( 16) is the linear conduction profile while the additional cosine term at O( 3 ) renders the thermal boundary conditions for θ 3 to be homogeneous (see Equation ( 14)).At O( ), the leading order disturbance equations are, These homogeneous equations form the linearised stability problem for the uniformly heated layer, and are to be solved subject to the following homogeneous Dirichlet boundary conditions, identical conditions will also apply for ψ i and θ i , i = 2, 3. We choose to use the following roll eigensolutions, where A = A(τ) is a complex amplitude which is a function of the slow timescale, τ.The complex forms used in Equation ( 21) allow for a straightforward modelling of phase variations in the solutions.
At O( 2 ) the equations for ψ 2 and θ 2 are, When the solutions given in Equations ( 21) are substituted into Equation ( 23) the latter becomes, and therefore ψ 2 and θ 2 are found to be, At O( 3 ) the equations for ψ 3 and θ 3 are, After the substitution of Equations ( 21) and ( 25), the above equations become, The inhomogeneous terms in Equations ( 28) and ( 29) contain components which are proportional to the eigensolution of the homogeneous form of the equations, namely, the solutions given in Equation (21).The value for R 2 is presently arbitrary but in general this system of equations cannot be solved unless R 2 takes a specific value.The derivation of a solvability condition will yield the desired amplitude equation for A in terms of R 2 , the satisfaction of which guarantees the solution of Equations ( 28) and (29) up to an arbitrary multiple of the eigensolution.
If Equations ( 28) and (29) are written in the form, and if we define the quantities, ψ * = iπe −iπx sin πz and θ * = 1 2 e −iπx sin πz, which are the coefficients of A in Equation ( 21), then the following identity may be formed: The left hand side of Equation (31) may be shown to be precisely equal to zero by using integration by parts, and thus the solvability condition is that, Upon application of this solvability condition to Equations ( 28) and ( 29) we obtain the following amplitude (i.e., Landau) equation, The coefficients of the last two terms of Equation ( 33) may be simplified by rescaling A, τ and R 2 in order to obtain the following canonical form of the amplitude equation, The scalings which are required to do this are, If we omit the asterisks for simplicity of presentation, then the complex amplitude equation governing weakly nonlinear convection is, where τ * has been replaced by t for simplicity of presentation from this point onwards, and this t should not be confused with that used in § §1 and 2. The rest of the paper is devoted to the properties of the real and complex solutions of this equation.For any given ω the forcing period is T = 2π/ω.

Initial Considerations and Context
Rees and Riley (1986) [12] obtained the following amplitude equation for a steady perfectly resonant imperfection, and this is equivalent to the setting of ω = 0 in the present paper.The steady solutions of Equation (37) are shown in Figure 2 and this Figure shows that no zero-amplitude solution is possible.When R 2 < 0 a unique solution is obtained, and it satisfies the asymptotic relation, A ∼ −R −1 2 , when R 2 has a large amplitude.But as R 2 increases, the amplitude of this solution branch also increases until it approaches asymptotically to A ∼ R 1/2 2 when R 2 is large and positive.The upper branch of this bifurcation diagram is stable.However, when R 2 > 3/2 2/3 , there exist two other solution branches, and both of these are unstable with respect to perturbations in phase when the porous layer is unbounded horizontally.The middle branch is also unstable unconditionally to perturbations in amplitude.37).The upper branch is stable.The other two branches are unstable with respect to perturbations in phase, while the middle branch is also unstable with respect to perturbations in amplitude.
An alternative way of visualizing the stability properties of Equation (37) which are described above is shown in Figure 3.In this Figure are displayed solution trajectories of A as a complex amplitude for different values of R 2 and for different initial conditions placed on a unit circle in the complex plane.For all choices of R 2 , the values of A converge onto the solution which is on the positive real axis.When R 2 is sufficiently large and the initial condition is close to the negative real axis, then the amplitude increases rapidly to close to R 1/2 2 and then the phase varies relatively slowly until the positive real axis is reached.The variation in phase is equivalent to a slow physical migration of the rolls towards the stable location, and this means that |A| has also been maximised.
Figure 3. Numerical solutions of Equation (37) for initial conditions placed on a unit circle in the complex plane, for Ra = −5, −1, 0, 1, 2 and 5. Showing how solutions evolve to the unique solution on the positive real axis.
The usefulness of the above picture lies in the fact that it gives a good insight into what might happen when ω is no longer zero.Thus when cos ωt > 0, Figure 3 gives an indication of the trajectories, while the mirror image about the vertical axis corresponds to when cos ωt < 0. Therefore we may draw an a priori conclusion that there might be the possibility of cases where there could be permanent oscillations in phase, given that a positive real value of A forms the attractive stable solution when cos ωt > 0 and its negative counterpart does so when cos ωt < 0.

Numerical solutions
We may expand Equation (36) by splitting it into its real and imaginary components: A(t) = B(t) + iC(t), where B(t) and C(t) are both real.This yields the following pair of coupled non-linear equations, and From Equation (39) we see that C = 0 is a possible solution, which corresponds to real values of A and, given the definition of A in Equations ( 21), convection is exactly in phase with the thermal forcing when both A and cos ωt are of the same sign.
In this subsection we shall concentrate solely on real solutions of Equation ( 38) with C = 0, namely, This equation was solved using the classical 4th order Runge Kutta method with at least 2000 timesteps within a period until all transients have decayed and a steady periodic state is achieved.Integration took place over a whole number of periods of the forcing, i.e., up to the time t max = N(2π/ω) where N is an integer.For each period, the solution was compared with that over the previous period, and when the maximum absolute difference over all the time steps was less than 10 −6 then convergence was deemed to have occurred.
In Figure 4 steady real periodic solutions are shown for a variety of frequencies, ω, and for a range of values of R 2 lying between −5 and 3.These curves are typical of all periodic real solutions.There are two different types of solution which occur as R 2 and ω are varied: (1) a single-signed oscillation, which we call Type I, (2) solutions taking both positive and negative amplitudes, which we call Type II.
It is clear from Figure 4 that, for all values of ω, single-signed solutions (Type I) always exist when R 2 is sufficiently large and that double-signed solutions (Type II) exist when R 2 is negative.However, the transition between these two states isn't straightforward to see in this Figure .When the forcing period is large (e.g., when T = 100) there appears to be a sudden transition between the two Types.However, when the forcing period is small (such as when T = 1 and T = 5) the solution curves simply descend as R 2 decreases.There will be a precise value of R 2 (which we will call R zero 2 ) at which the minimum value of B is zero; these values are listed in Table 1.But the consequence is that there appears to be an intermediate range of R 2 which seems to be transitional between the two Types.
Further analysis of these curves indicates that the characterisation of the two Types based on whether a particular curve is single-signed or double-signed is too simplistic.A much better alternative is to consider the variation in the maximum and minimum values of B over one period; such curves are displayed in Figure 5. Concentrating first on the continuous lines for T = 5 and for low values of R 2 , we see that B max + B min = 0, and therefore the corresponding solution curves for B are double-signed and are of Type II.Once R 2 reaches roughly 0.838 there is a clear bifurcation to a new form of solution, and it is this which should be labelled as being of Type I. Again concentrating on the continuous curves, we see that B max and B min still have the opposite sign for a further small range of values of R 2 , but once R 2 exceeds roughly 1.043 the solution then becomes single-signed.
In Figure 5 the continuous curves correspond to the use of B(t = 0) > 0 as the initial condition, whereas the dotted lines correspond to B(t = 0) < 0, or, equivalently, to solutions with a phase-shift of π compared with the solutions with the positive initial condition.Inclusion of these latter curves give rise to sets of curves which are more obviously of supercritical bifurcation type, and it is this point of bifurcation which marks the exact transition between solutions of the two Types.Although we shall not demonstrate it, Type II solutions also exist above the bifurcation point and have been computed, but they are unstable.Given that the bifurcation point is where the Type II solutions become unstable, we may determine these points numerically.If we set A = B + D into Equation (40), where D is both real and asymptotically small, then we obtain the following linearised equation for D: The bifurcation point now corresponds to nonzero solutions for D for which D(0) = D(2πω).This was incorporated into the above-mentioned shooting method code where, for a chosen value of ω, the critical value of R 2 may then be found.Although a Newton-Raphson scheme was used, solutions converge only linearly, and with increasing difficulty as ω becomes small.However, these critical values of R 2 , named R bif 2 are also shown in Table 1.Now we proceed to the analysis of the high frequency cases, and Figure 4 shows clearly that the solutions are of Type I and are essentially constant but with a small periodic variation superimposed.It proves convenient to use a new time scale, τ = ωt, in Equation (40); this yields the equation, Equation ( 42), when expanded using the following series, gives a sequence of equations for the real functions, B 0 (τ), B 1 (τ) and so on.Given the above observation regarding the form of the solution, we will expect that B 0 is constant.
At O(ω) in the expansion we have the equation, B 0τ = 0, which confirms the numerical evidence that B 0 is a constant.Therefore at O(1) in the expansion we have the equation, This is one equation for two unknowns but, given that B 0 must be constant, the solutions for B 0 and B 1 are, No other solutions are possible which are physically reasonable, for if B 0 were any constant other than √ R 2 (or alternatively − √ R 2 ) then B 1 would contain a term which grows in time, and this possibility must be discounted given the evidence of Figure 4.The B 1 term shows that the leading unsteady response is π/2 out of phase with the cosine forcing when the forcing frequency is high.
At O(ω −1 ) the governing equation is and upon substitution of the previous solutions we have Thus, the solution for this equation is, where α is a constant which is determined later.At O(ω −2 ) the equation is (1) and these are described in detail below.In the following analysis we shall continue to use Equation (42) for mathematical convenience.We begin by expanding the solution of Equation (42) using the following power series in ω, At O(1), the resulting equation for B 0 is, As there is no time-derivative in this equation B 0 varies smoothly and quasi-statically when either R 2 < 0 or R 2 > 3/2 2/3 ≈ 1.889881.That this is so may be seen by considering the quasi-static variation of B 0 as depicted in Figure 6 where the various symbols represent the solution for three representative values of R 2 .The curves displayed in each frame of Figure 6 represent all the possible solutions of (57) at different values of τ over half a forcing period.In each frame, and from left to right, the symbols indicate potential solutions for (i) R 2 < 0, (ii) 0 < R 2 < 3/2 2/3 and (iii) R 2 > 3/2 2/3 .The solid circles depict how one might predict how an unsteady solution might evolve from the initial conditions given when ωt = 0, while the circles depict possible solutions but these do not arise with the chosen initial conditions.
When R 2 < 0 there is only one possible solution branch and therefore B 0 can vary quasi-statically as indicated by following the left hand circles in each subfigure.When R 2 > 1.889881, which is represented by the right hand symbols in each frame, there are always three branches and it is again possible for each outer branch solution to vary quasi-statically and remains on that branch.In this case we observe that the upper branch always exists and therefore there is no mechanism to cause the solution to change to another branch.

Sudden Transitions Between Branches
However, when 0 < R 2 < 1.889881, then there are parts of the forcing period when three branches exist, and parts when only one exists.If we follow the movement of the filled circle from frame 1 (ωt = 0) onwards in Figure 6, then we see a gradual decrease in the amplitude of the solution on the upper branch until frame 4 (ωt = 3 4 π).At a point in time somewhere between this and ωt = π, the upper branch ceases to exist at that value of R 2 , and the only branch which exists is then the lower branch.Therefore we are forced to expect to see a rapid change in the solution when τ is close to those values corresponding to when R 2 is at the nose of the conjoined branches.Examples of this have already been seen for T = 100 and R 2 < 1.4 in Figure 4.It is straightforward to show that this 'quasi-static' transition happens when where cos τ * = −2(R 2 /3) 3/2 and for integer values of n.One example of the sensitivity of the solution curves to the value of R 2 is shown in Figure 7 where ω = π/100 (T = 200) has been chosen as the representative small value of ω.At this value of ω the transition between flows of the two types occurs between R 2 = 1.871665 and R 2 = 1.871666, and therefore there is an obviously sharp transition between solutions of Type 1 and of Type II.Solutions for the neighbouring values, R 2 = 1.4,1.8 and 2 are shown in Figure 8, together with the corresponding quasi-static solution curves, for comparison.For the Type II solutions depicted, the true solution is given very accurately by the quasi-static solution, but as it approaches the nose of the quasi-static curve, the value of B drops precipitously to the lower branch and is then given very accurately by the negative quasi-static solution.We now proceed to a study of the process of transition from one branch to another.On letting τ = τ * + ω t in Equation ( 42), where τ * is defined in Equation ( 59), we have Thus t is simply a phase-shifted version of t.The function cos(τ * + ω t) may be expanded for small values of ω according to, where S = cos τ * = −2(R 2 /3) 3/2 < 0. Thus, if Equation ( 60) is rewritten in terms of the leading terms given in Equation ( 61), then we obtain, Upon rewriting R 2 in Equation (62) in terms of S, and rescaling using then we obtain, The right hand side of Equation (64) may be simplified, and if also rename the last term to be Sω t then the equation becomes, where S = 2(2/3) 5/2 ( √ 1 − S 2 )/S, which also depends on R 2 .We may now use the method of dominant balance to determine suitable scalings for further analysis.By setting, B = 1 + ω a B and t = ω b t.It was found that a = 1/3 and b = −1/3 form the suitable powers of ω for the subsequent analysis.Hence, the substitution of these values into Equation (65) gives the following equation, Figure 9 shows the numerical solution of Equation (66) the case when S = −0.5 where the solution separates from the quasi-static solution before taking a sharp descent.The shape of the tranient solution shown in Figure 9 has obvious qualitative similarity with the Type II solutions shown in Figure 8.This asymptotic solution is also fully representative of all other values of S since S may also be scaled out of Equation ( 66) using yet another transformation.The solution shown in Figure 9 becomes negatively infinite at a finite time.Given that the present analysis is a local analysis, we conclude that this is equivalent to jumping from the upper branch to the lower branch, which is in accord with many of the curves for small values of ω (or large values of T ) which are shown in Figure 4. Upon going back through the various transformations, this indicates that the transition timescale is of order ω −1/3 in duration, which is very short compared with the forcing period, 2π/ω.
An identical analysis may be used to describe the reverse transition from a lower branch to an upper branch when cos τ * = 2(R 2 /3) 3/2 and 3 2 π < τ * < 2π, i.e., exactly half a period after the first transition; see Equation (59).

Stability Analysis for Real Solutions
So far we have concentrated on analysing solutions for which A is real, where the convective pattern remains stationary relative the spatial pattern of the temperature variations on the boundary of the layer.The work of Banu and Rees [22] (2001), which studied a moving resonant wavelength thermal forcing, naturally found cellular patterns which also move along the layer.Mathematicially such movement is accounted for by allowing the amplitude, A, to be complex.
In the present context it is possible to predict the existence of complex solutions by considering again the sketches in Figure 6 If we look at the R 2 = 4 solution (filled circle) which is on the upper branch, then once τ is as large as that represented in frame 4 (ωt = 3π/4), that solution is unstable with respect to perturbations in phase.There is therefore an a priori expectation that the cellular will move towards the closer position where the heat transfer will be maximised, Therefore it is to be expected that the phase of the pattern should be a function of time, in general.

General Linear Stability Analysis
We may determine the linear stability equations for the real solutions, B(t), given above, by substituting into Equation (36) and by linearising, given that we are setting |D| 1.This gives the following equation for the small-amplitude disturbance, D: We may split D into its real and imaginary components using D = D R + iD I , and therefore we have the following equations, and The form of these decoupled equations makes it clear that imaginary disturbances will always become unstable at a smaller value of R 2 than real disturbances will, and therefore our linear theory will consider only imaginary disturbances.We also note that Equation (77) is identical to Equation (41), and therefore we may already conclude that real solutions of Type I, which exist after the supercritical bifurcation, will always be unstable to phase disturbances.

Numerical Solutions
The basic state which we analyse for stability is given by the solution of Equation (40) while the linearised disturbances satisfy Equation (78) above.There are various ways in which stability criteria may be found numerically.One of these involves integrating both Equations ( 40) and (78) forward in time using a non-zero initial condition for D I .If we denote by D n the value of D I after n forcing periods, then the critical value of R 2 corresponds to when the following relation is satisfied, One alternative procedure would be to use a shooting method based on the classical fourth order Runge-Kutta code mentioned earlier where the boundary condition, D I (0) = D I (2π/ω), is imposed in order to find the critical value of R 2 .A second alternative relies on the fact that D I may be found analytically in terms of B: If there is no overall growth over one forcing period then neutral conditions exist and, in terms of τ, the following condition needs to be met: this integral may also be encoded within a shooting method code.
The results of an extensive set of calculations are shown in Figure 12 and Table 2 where the numbers given are correct to 8 decimal places.The detailed numerical values indicate that When R 2 and ω take values which are below the stability curve, which is the dotted curve in Figure 12, then real solutions are stable.Otherwise, the imaginary part grows and the final periodic state is complex, implying that the convection cells move horizontally.This neutral curve also lies below the whole of the curve marking the transition from Type II solutions to Type I solutions as R 2 increases, as mentioned above.This means that real type I solutions are unstable when the porous layer is of infinite horizontal extent.On the other hand, if the porous layer were to be of finite width, and specifically a whole number of spatial periods of the boundary forcing, then disturbances which alter the phase of the disturbance uniformly, as we have used here, will not satisfy any acceptable sidewall boundary conditions, and therefore real solutions which are stationary will remain stable, at least within the weakly nonlinear regime.

Stability for Large Frequencies (ω 1)
In the large-ω limit it is possible to supplement the above numerical data using the asymptotic basic state solutions given in Equations ( 53)-(55).For the upper branch solution given in Equation (53), substitution into Equation (81) yields, which shows that this solution is always unstable.An alternative way of viewing Equation ( 83) is that the imaginary disturbance takes the form of a time-periodic function multiplied by e τ/ω 2 , and therefore the disturbance grows very slowly in time.An almost identical analysis shows that the lower branch solution given in Equation ( 54) is also unstable unconditionally.
If we take the leading term in the middle branch solution Equation (55), then we have This solution will be unstable once R 2 exceeds 1/2ω 2 , and therefore the one-term large-ω expression for the neutral curve is Although this neutral curve compares quite well with our numerical data, it is worth attempting to find the next term in the large-ω expansion.Given that powers of R 2 appear in Equation (55), where R 2 is assumed to be of O(1), it will be necessary first to rescale R 2 in order to carry out a rigorous large-ω analysis.Therefore we set, and we are now solving,

Complex Solutions
Equations ( 38) and (39) are the expanded forms of Equation (36) which is both complex and non-linear.Based on the stability results given in Table 2 of the linear form of Equations ( 38) and (39), we could begin to seek the qualitative behaviour of solutions at higher Rayleigh numbers.

Numerical Simulations
A Runge Kutta code was run until a steady periodic complex solution was achieved.The solutions thus obtained are independent of initial conditions and transients decay quite quickly when the real initial condition is taken to be √ R 2 .A nonzero imaginary part is used to seed complex solutions, for if a zero imaginary part were to be used then the solution will always remain real.When ω is small the steady periodic state is attained after 2 or 3 forcing periods.On the other hand, when ω is large, it takes many periods to achieve a steady periodic state; this is consistent with the very small exponential growth rate which is a consequence of Equation (83).A set of different R 2 values for each selected values of ω is used in the code.A selection of these solutions are shown in Figure 14, and all represent a distinctive pattern showing oscillations above the real axis whenever R 2 exceeds its critical value.
An interesting observation is found that the stable solutions oscillate about the zeros on the real axis while the unstable real solutions become closed trajectories on the complex plane.For the case of unstable solutions, when cos ωt is positive the steady solution is attracted to the stable branch therefore, as the result the solution curve moves closer to the real axis where it follows the quasi-static solution.This result is especially obvious when the value for ω is small.When convection changes sign or rather when the boundary forcing changes sign, which is analogous to cos ωt being negative, the stable branch changes position after going through the supercritical pitchfork bifurcation.Hence, the solution curve jumps to the second branch where it is stable before it follows the quasi-static solution again and complete the phase.The solutions are therefore unstable to variations in phase which is clearly seen in Figure 14 when ω values are small.

Solutions for Large Values of R 2
Figure 14 shows that solutions for large values of R 2 tend to show quite small oscillations about a value close to A = √ R 2 i.Therefore, we shall conclude our analysis of Equation (36) by considering the large-R 2 limit of its solutions.The analysis begins by adopting the following series expansion, On substituting into Equation (36), the equation which arises at In general, any complex value of A 0 satisfying |A 0 | = 1 is a solution of this equation, but the fully numerical solution suggests that A 0 = i should be set.
At O(R 2 ), the equation we obtain is, The only valid solution of Equation ( 95) is a real function of t, but it must allow for the possibility of periodic oscillations in the real part of A. However, this equation is silent on what that function is and therefore we shall set A 1 = a 1 (t) for now and determine a 1 later.
At O( √ R 2 ) we have, Therefore the solution for A 2 must be a combination of real and imaginary parts, and we obtain, where a 2 (t) is also an unknown real function at present.At O(1) we obtain, which, after substitutions and simplification, reduces to the form, The substitution of A 3 = a 3 (t) + ib 3 (t) into this latest equation above reduces it to Since a 1 is real, the real and the imaginary parts may be separated to give, The real part gives us, which provides the required oscillation in the A 1 term.Therefore we may now write, where at present both a 2 (t) and a 3 (t) remain unknown.
At O(1/ √ R 2 ), the equation to solve is, and, with simplification, this becomes, Upon using the substitution given in Equation (97), then the real part of Equation (105) becomes, The general solution of this equation is where a 3 (t) and a 4 (t) are now the remaining unknowns.Further terms have been obtained in the same manner although the algebra becomes very lengthy and therefore it is omitted.The expression for a 3 is found at higher order in the expansion and it turns out that a 4 = 0. To summarise, we have obtained the following asymptotic solution for Equation (36), Comparisons of the asymptotic solution and the numerical solution are shown in Figure 15.The asymptotic solutions show an increasingly good quantitative agreement with the numerical solution as R 2 increases.Clearly the large-R 2 state consists of cells which are π/2 out of phase with the boundary forcing, and which twitch slightly in the horizontal direction.0.6 0.4 0.2 0.2 0.4 0.6 1.9

Conclusions
We have studied the behaviour of two-dimensional convective rolls in a horizontal fluid-saturated porous layer heated from below, taking into account the effect of time-periodic resonant thermal boundary imperfections of small amplitude.A weakly non-linear theory has been used to derive the governing amplitude equation and to analyse the stability of the motion.
There exist different solutions of this amplitude equation in different regions of Rayleigh number/wavenumber space.If the flow is constrained to be stationary, which will happen if the layer is bounded by suitably-placed sidewalls with impermeable and insulating boundary conditions, then what is termed a Type II flow exists and is stable when Ra is below a frequency-dependent critical value.As Ra increases there is a supercritical bifurcation to what is termed a Type I flow.Flows of Type II are double-signed, whereas Type I flows are usually single-signed-exceptions are situated close to the bifurcation point.Detailed asymptotic solutions are presented for small frequencies where the numerical solutions display localised and quite rapid changes in sign.
When the layer is unbounded horizontally, then the real-A stationary solutions are susceptible to disturbances in phase, and we find that all Type I solutions are unstable to this type of disturbance.Type II real solutions may also be destabilised if R 2 is sufficiently large to become one with complex amplitude, which corresponds to a convection pattern which oscillates horizontally.
Numerical solutions have in many cases been supplemented by asymptotic solutions and these have shown excellent agreement.
Three dimensional disturbances have not been considered in this paper.Given that the boundary disturbances have a spatial pattern which has the critical wavenumber, it is our a priori expectation that three-dimensional effects are subdominant.However, we think that three-dimensional effects will be important when the spatial wavenumber of the disturbance is different from k = π.Certainly it is the case that such patterns are dominant when these boundary forcing terms are steady in time; see Rees and Riley 1989(a,b) [12].
Finally we note that some of the work described here on moving cells may need to be modified if a layer of long but finite length were to be considered.It has already been stated that disturbances in phase cannot arise in a short layer due to the presence of sidewalls.However, if the layer is of finite but sufficiently large extent (which will be of O(δ −1/3 )), then it may well be possible for cells to be compressed slightly near those sidewalls, and this would allow for phase movements in the bulk which then reduce in magnitude as the sidewalls are approached.In these the resulting time-dependent flow may be visualized by having stationary cells fixed or pinned to the two sidewalls, rather than being generated or annihilated, and then the magnitude of the time-variation of the phase of all other cells will increase to a maximum midwall between the sidewalls.If the layer is sufficiently long, then there will be a substantial region in the middle of the layer where the present theory is valid.
It is also intended to extend the present work well into the nonlinear regime.

Figure 1 .
Figure 1.Neutral curve for the classical Darcy-Bénard problem with n = 1; see Equation (1). the critical values are given by Ra c = 4π 2 and k c = π.

3 AFigure 2 .
Figure 2. The steady solution curves for the amplitude equation Equation (37).The upper branch is stable.The other two branches are unstable with respect to perturbations in phase, while the middle branch is also unstable with respect to perturbations in amplitude.

Figure 4 .
Figure 4.The effect of different Rayleigh numbers on real solutions of Equation (36) for the following selection of forcing periods: T = 1, 5, 20 and 100.The respective frequencies are ω = 2π, 0.4π, 0.1π and 0.02π.The uppermost curve corresponds to R 2 = 3; the dash-dotted curve to Ra = 0 with intermediate curves corresponding to intervals of 0.2 in Ra; the dotted curves correspond to Ra = −1, −2, −5 and −10.

Figure 5 .
Figure5.The bifurcation diagram for the given values frequencies/periods.The curves correspond to B max and B min over period.The dashed lines correspond to solutions for which B(t = 0) is negative, while continuous lines have B(t = 0) > 0.

Figure 6 .
Figure 6.Depiction of the bifurcation diagram corresponding to quasi-static real solutions of Equation (40) over half a period.The five frames correspond to ωt = 0, π/4, π/2, 3π/4 and π.In each frame, and from left to right, the symbols indicate potential solutions for (i) R 2 < 0, (ii) 0 < R 2 < 3/2 2/3 and (iii) R 2 > 3/2 2/3 .The solid circles depict how one might predict how an unsteady solution might evolve from the initial conditions given when ωt = 0, while the circles depict possible solutions but these do not arise with the chosen initial conditions.

Figure 7 .
Figure 7.An example of a catastrophic change in the qualitative nature of the solution for the small frequency, ω = π/100.The solid line corresponds to R 2 = 1.871666 and the dashed line to R 2 = 1.871665.

Figure 12 .
Figure 12.Showing the value of R 2 as a function of ω at which solutions of Type II undergo a supercritical bifurcation to solutions of Type I (continuous line), and the neutral stability curve with respect to imaginary disturbances (dashed line).

4 Cω 5 Figure 14 .
Figure 14.The effect of different values of R 2 and ω on solution trajectories.Bold values of R 2 indicate values for which solutions are real.All solution trajectories follow that of the arrow shown in the top right subfigure.

Figure 15 .
Figure 15.Continuous curves show the solution trajectories for the given values of R 2 .The left hand column corresponds to ω = π/2 and the right hand column corresponds to ω = π/20.Dashed lines indicate the large-R 2 asymptotic solution.

Table 1 .
The values of R 2 corresponding to the transition between Type I and Type II solutions of Equation (40).R zero 2 marks where B min = 0, while R bif 2 marks the bifurcation point in Figure5.

Table 2 .
Numerical solutions showing the relationship between ω and R 2c corresponding to neutral stability with respect to phase perturbations.In the right hand subtable we also show R 2c,asymp from Equation (92) and its error.
)which decays exponentially in time.Therefore we may set a 2 = 0.The solution proceeds by substituting A 4 = a 4 (t) + ib 4 (t) into Equation (105) and this gives,