Linear and nonlinear stability of Rayleigh–Bénard convection with zero-mean modulated heat ﬂux

Linear and nonlinear stability analyses are performed to determine critical Rayleigh numbers ( Ra cr ) for a Rayleigh–Bénard convection conﬁguration with an imposed bottom boundary heat ﬂux that varies harmonically in time with zero mean. The Ra cr value depends on the non-dimensional frequency ω of the boundary heat-ﬂux modulation. Floquet theory is used to ﬁnd Ra cr for linear stability, and the energy method is used to ﬁnd Ra cr for two different types of nonlinear stability: strong and asymptotic. The most unstable linear mode alternates between synchronous and subharmonic frequencies at low ω , with only the latter at large ω . For a given frequency, the linear stability Ra cr is generally higher than the nonlinear stability Ra cr , as expected. For large ω , Ra cr ω − 2 approaches an O ( 10 ) constant for linear stability but zero for nonlinear stability. Hence the domain for subcritical instability becomes increasingly large with increasing ω . The same conclusion is reached for decreasing Prandtl number. Changing temperature and/or velocity boundary conditions at the modulated or non-modulated plate leads to the same conclusions. These stability results are conﬁrmed by selected direct numerical simulations of the initial value problem.


Introduction
The convection configuration considered in this paper has emerged from the study of radiatively driven convection in ice-free Lake Superior.During springtime warming of this freshwater lake, observational data in Austin (2019) appear to show that an instability arises each day near the surface and carries heat through the entire water column on a time scale of hours.The temperature of the lake is between approximately 0 • C and 4 • C, which means that the water is in the anomalous regime where increasing temperature increases rather than decreases density.
The instability can be understood physically as follows.Observations show that the water column begins the day at a uniform temperature throughout.As the sun comes up, radiative heating penetrates into the water column, with the heating concentrated near the surface and falling off exponentially with depth.Because the water is in the anomalous regime where temperature increase leads to density increase, heating at the surface increases the density of the water there.Buoyancy then causes the denser water to sink.If the buoyancy effects outweigh the restraining effects of heat diffusion and viscosity, then an instability arises.In many bodies of water, radiative heating penetrates into only a small fraction of the water column, and here we therefore treat the limit of radiative heating confined to an infinitesimally small layer near the surface, meaning that we specify a time-varying heat flux at one of the boundaries rather than including a radiative source term in the governing equations.The two infinite surfaces with heat flux imposed at one boundary and temperature imposed at the other makes this essentially a Rayleigh-Bénard (RB) configuration, but with an imposed flux that is modulated in time rather than an imposed steady temperature difference.
Most previous work has considered modulation on top of a background temperature gradient.Here, we treat modulation with zero mean, meaning that the average of the temperature difference between the top and bottom surfaces is zero for the unperturbed base state.With zero-mean modulation, if the amplitude of the boundary heat-flux modulation is set to zero, then nothing interesting happens as the water column is stably stratified from gravity and uniform in temperature.Only with a non-zero modulation amplitude is there a possibility for an unstable configuration.
The linear stability of modulated convection has been studied in earnest since at least Gershuni & Zhukhovitskii (1963), who looked at the low-frequency limit of modulated temperature gradient in standard RB convection.Just as temperature boundary conditions were considered first in standard RB convection, temperature modulation was considered first for modulated RB convection.In particular, the combination of no-stress velocity conditions and imposed temperature boundary conditions allows an analytical solution in terms of sine functions to be obtained, which was the approach taken in Venezian (1969), where the amplitude of modulation of the boundaries in standard RB convection was taken as a small parameter.
These were followed by a number of studies on the linear stability of modulated convection.Few authors, however, addressed zero-mean modulation, with the exception of Yih & Li (1972), Gershuni & Zhukhovitskii (1976), Or & Kelly (1999) and Souhar & Aniss (2016).These authors investigated boundary temperature modulation, but no one appears to have addressed boundary heat-flux modulation.Davis (1976) reviewed the stability (linear and nonlinear) of a variety of time-periodic flows, including thermal instabilities, but did not explicitly mention zero-mean modulated flow or heat-flux modulation.
In addition to linear stability, we also examine nonlinear stability of modulated convection using the energy method.The first major work using the energy method to establish nonlinear stability in fluid dynamics appears to be Joseph (1976), though it was used before that by a variety of researchers, including Serrin (1959), who cited Reynolds and Orr as the progenitors.More recent textbooks covering the energy method include Doering & Gibbon (1995) and Straughan (2004).The work most relevant to our concerns is Homsy (1974), which investigated gravity modulation as well as modulation of the boundary temperatures.
To summarize the two approaches to determining stability, linear stability analysis establishes a sufficient condition for instability, in this case a Rayleigh number above which at least one infinitesimal disturbance grows exponentially in time.Nonlinear stability analysis establishes a sufficient condition for stability, in this case a Rayleigh number below which the energy of all disturbances eventually decays.In the present case with a time-periodic base state, we consider two possibilities.For asymptotic stability, the disturbance may grow during a cycle but overall experiences net decay, whereas for strong global stability the disturbance decays exponentially in time.
In the present paper, we consider convection in a layer of fluid that is infinite in the horizontal directions.We investigate zero-mean modulation of the heat flux at the bottom boundary of a standard fluid layer, which, from the symmetry of the modulation profile, is equivalent to modulation at the top of a water layer in the anomalous 0-4 • C regime, as would be the case for a lake.For comparison, we also give results for zero-mean modulation of the temperature at the bottom boundary, though we do not detail the derivation of the equations.After discussing the set-up and governing equations in § 2, we go through the calculation of linear and nonlinear stability thresholds in § § 3 and 4, respectively.We then present results in § 5. Our main conclusions are further confirmed by selected direct numerical simulations (DNS) of the initial value problem, presented in § 6.We finally discuss our results and future works in § 7.

Set-up
We consider two parallel plates extending infinitely far in the horizontal x-and y-directions, containing fluid satisfying the Boussinesq equations, where asterisks represent dimensional quantities.Here, u * is the velocity, T * is the temperature measured with respect to the reference temperature at the upper boundary, ρ 0 is the density at the reference temperature, g is the acceleration due to gravity, ν is the kinematic viscosity, κ is the thermal diffusivity, α is the thermal expansion coefficient, and p * is the pressure.Figure 1 shows a schematic of the configuration.The velocity boundary conditions can be either no-stress or no-slip, and the temperature boundary conditions are (2.4) where d is the domain size in the vertical z-direction pointing up, H is the amplitude of the modulated heat flux, and k is the conductivity.We non-dimensionalize using and an appropriate scaling for pressure.We note that choosing to scale time by a diffusive time scale (e.g.d 2 /κ) may be more appropriate in certain cases such as in the low-frequency limit ω * → 0.
Using the non-dimensional variables in (2.6a-d), the non-dimensional governing equations become (2.9) The temperature boundary conditions become (2.11) The non-dimensional frequency ω, Rayleigh number Ra, and Prandtl number Pr, are defined as We write the base state velocity, temperature and pressure as u B , T B and p B .For the stability analysis, we take a base state with no motion (u B = 0) and a temperature profile satisfying (2.9) with the velocity set to zero, namely (2.13) and the boundary conditions (2.10) and (2.11).We write the solution for the base state temperature as

Linear stability calculation
For linear stability, we consider small perturbations to the base state and write Using these in the governing equations (2.7) and (2.8), neglecting products of perturbations, and isolating the vertical velocity component (e z • u p = w p ) results in where The governing equations have constant coefficients in space, therefore the horizontal spatial components of the functions can be analysed using normal modes, so that we can write the perturbations as w p (x, t) = e ik x x e ik y y w(z, t), (3.4) θ p (x, t) = e ik x x e ik y y θ(z, t). (3.5) The resulting equations are where To manage the z dependence, we use Chebyshev polynomials.One possibility is to use Chebyshev differentiation matrices, as discussed in Weideman & Reddy (2000) and Trefethen (2000), so that w and θ are solved for at specific grid points.Another possibility is to express w and θ as Chebyshev polynomials and then use collocation or Galerkin projection to remove the z dependence so that the coefficients in the Chebyshev expansion of w and θ become the relevant unknowns, which is the approach taken in Or & Kelly (1999).The boundary conditions must be satisfied in each case.For details on the numerical methods used, see Appendix A.
The resulting matrix equation can be written as with x = (w, θ) T , where w and θ are vectors of coefficients or of the solutions at chosen grid points, depending on the method chosen.The base state gradient is expressed as ∂ z T B = T 1 (z) e it + T −1 (z) e −it and discretized using the same method as the solutions.
The system of ordinary differential equations (ODEs) represented by (3.8) has coefficients that are periodic in time, and we therefore use Floquet theory.There are two ways in which we can use Floquet theory: the Floquet-Fourier-Hill (FFH) and monodromy matrix methods.The FFH method requires solving an eigenvalue problem, while the monodromy matrix method requires solving a system of ODEs.Here, we use the FFH method because it is more efficient computationally.For details on the FFH method, see Deconinck & Kutz (2006).

Floquet-Fourier-Hill method
For the FFH method, we use Floquet theory to decompose w(t) and θ(t) into an exponential function of time multiplying a function that is periodic in time with the same period as the coefficients (2π).The solution vector is then x n e int , (3.9) where μ is the Floquet exponent.Using this representation in (3.8) leads to Orthogonality of the exponential functions leads to (3.11) which is an infinite system of coupled equations for the Fourier coefficients and can be treated as an eigenvalue problem for Ra or μ.We solve this coupled set of equations numerically by truncating the Fourier series and solving the resulting generalized eigenvalue problem, which is block tridiagonal on one side and block diagonal on the other.
The eigenvalue problem depends on Pr, ω, k, Ra, μ, the number of Fourier modes, the number of Chebyshev modes or grid points, and the boundary conditions.Once other parameters are fixed, either the Rayleigh number Ra or the Floquet exponent μ can be considered as the eigenvalue.The critical Rayleigh number for given Pr and ω is the lowest Rayleigh number found through varying k that results in Re(μ) = 0.
If the Rayleigh number is treated as the eigenvalue, then we fix Re(μ) = 0 to look for marginal stability.If we write μ = μ 0 + im in (3.9), with 0 Im(μ 0 ) 1 (where Im(•) indicates 'imaginary part of') and m a positive integer, then we find x n e i(n+m)t . (3.12) We see that m serves only to shift the association between the coefficient index and the frequency of the exponential that it accompanies, and we can therefore restrict Im(μ) to lie between 0 and 1 without loss of generality.The lowest Rayleigh number found over all wavenumbers is the critical Rayleigh number.If, instead, the Floquet exponent μ is chosen as the eigenvalue, then k and Ra must be swept through to find the minimum value of Ra resulting in Re(μ) = 0.
For our numerical results, we have generally found the critical Rayleigh number by treating Ra as the eigenvalue.We have then checked the resulting critical Rayleigh number and wavenumber by using these values and treating μ as the eigenvalue to ensure that Re(μ) is truly close enough to zero to represent the stability threshold.Furthermore, we have checked the surrounding (k, Ra) parameter space to be sure that the critical Rayleigh number found is truly a local minimum leading to marginal stability.
We have generally used 18 Chebyshev grid points.The highest Fourier mode used when solving for Ra as the eigenvalue varied between 15 and 30, depending on the frequency, with more Fourier modes being necessary to reach a converged solution at lower frequencies.When solving for μ as the eigenvalue, we have been able to use sparse eigenvalue routines that use the fact that |μ| is small, which has allowed us to use a largest Fourier mode of 35 to 50.This is not possible for Ra values that are not small.

Low-frequency limit
At a certain O(1) value of ω, the eigenvalue problem resulting from the FFH method becomes too ill-conditioned to continue working.For small ω, the leading order of the z-derivative of the base state in (2.14) becomes independent of space and can be written as (3.13) This coincides with what we expect physically since for very slow modulation, the temperature profile will be almost linear, and its slope will be that imposed at the boundary, namely cos t.With the spatial dependence eliminated, we can compare the eigenvalue problems for stability in the modulated case to the eigenvalue problems for stability in the steady case with boundaries held at different temperatures, which has non-dimensional temperature profile slope −1.
For temperature modulation with no-stress boundary conditions, the governing equations reduce to a Mathieu equation in this limit, and a WKB analysis can be done, as discussed in Or (2001).For all other cases, a similar approach leads to a system of coupled Mathieu equations, which can be solved formally with an extension of the WKB ansatz to systems of equations.Unfortunately, connecting the WKB solutions through turning points is much more difficult with higher-order systems of equations than it is for the standard second-order ODE.Mathematical details are discussed in Wasow (1985), but there does not appear to be a simple way to use the WKB approach for the cases considered here.

Nonlinear stability calculation
To determine the threshold for nonlinear stability, we use the energy method and two notions of nonlinear stability: strong global stability and asymptotic stability, as used and described in Homsy (1974) for non-zero-mean temperature modulation.The analysis in Homsy (1974) uses the mean of the boundary temperature difference as the temperature scale, which is not possible for zero-mean modulation, and the base state is different, but otherwise the approach is similar.We therefore give only the essentials and refer the reader to Homsy (1974), Joseph (1976) or Straughan (2004) for more details.
The first step is to form an energy functional using power integrals.The first integral comes from taking the dot product of (2.7) with u and integrating over the volume.The second integral comes from writing the temperature as the base state plus a fluctuation of arbitrary size, T(x, t) = T B (z, t) + θ(x, t), using this expression in (2.9), and then multiplying by θ and integrating over the volume.Finally, we multiply the temperature integral by λ Ra and add the result to the momentum integral, where λ is a coupling parameter that we can later tune to achieve better stability results.The resulting equation for the time evolution of the energy is where R = √ Ra, φ = θ √ λ Ra, and the norms are We now define (4.3a-c) so that we have From (4.4), we develop strong global stability and asymptotic stability.

Strong global stability
For strong global stability, we divide both sides of (4.4) by D to find where H is the space of divergence-free functions satisfying the boundary conditions.We define 1 ρ λ (t) ≡ max where ρ λ (t) is periodic with the same period as the base state temperature gradient.We then define to arrive at From Poincaré's inequality, we have D α 1 E, with α 1 0. We therefore obtain For R < R S,λ , the energy decays exponentially or faster in time, which we call strong global stability.
To find R S,λ , we must first solve the variational problem for ρ λ (t) in (4.6), which upon using variational calculus with a Lagrange multiplier for the incompressibility constraint (see Straughan (2004) and Christopher (2021) for details), leads to We use normal modes and write The equations then become This is a generalized eigenvalue problem for ρ λ (t), which can be written as We solve this generalized eigenvalue problem by discretizing on a Chebyshev basis (as discussed in the linear stability section).To find R S,λ , we then minimize ρ λ (t) over time as specified in (4.7).Finally, we vary λ to find the best stability result, namely the highest R S,λ , which we define as the threshold for strong global stability, R S .Altogether, this is The eigenvalue problem for the linear stability threshold in standard RB convection can be written as The spatial operators in the two cases are equivalent, so ρ λ (t) must satisfy where each R j,steady satisfies the eigenvalue problem in (4.19) and (4.20).Rearranging, we have for each possible j.Now we use normal modes and perform R S = max λ min k min t ρ λ (t; k).Minimizing in time clearly means that we must take cos t = 1.We are then left to maximize over λ, which leads to λ = 1.Finally, we minimize the resulting ρ λ (t; k) over k, which leads to Because R L,steady ≡ min k {R j,steady }, we conclude that Ra S → Ra L,steady as ω → 0. This is the limit approached by the numerical results, as seen in figures 2 and 5 for modulated flux and temperature, respectively.

Asymptotic stability
For asymptotic stability, we start from (4.4) and write First, we define where H is the space of divergence-free functions satisfying the boundary conditions.This leads to To find ν λ (t), we use the variational calculus with a Lagrange multiplier for incompressibility, leading to The eigenvalue ν λ (t) is periodic, and we define the configuration as asymptotically stable if the integral of ν λ over one period is less than zero, defining We use normal modes as in (4.12) and find the generalized eigenvalue problem (4.31)We then discretize in z using a Chebyshev basis, and solve this generalized eigenvalue problem numerically in order to estimate νλ , the integral of ν λ (t) over one period.For fixed Ra, Pr and ω, we first sweep through wavenumbers and take the worst-case (largest) value for ν λ (t) at the chosen points in time over one period, which upon integrating in time gives us νλ .For efficiency, we have used Gauss quadrature with 30 grid points for the integral.Checks using more advanced integration methods indicate a relative error in the resulting Rayleigh number of well below 1 % when using Gauss quadrature with 30 grid points.We then vary λ to minimize this integral, and we define the result as (4.32) Finally, we find the largest R satisfying ν < 0 and define it as R A , so that

Low-frequency limit
For heat-flux modulation, we have not found any simplification to be possible in the low-frequency limit for asymptotic stability, and we therefore discuss temperature modulation only for symmetric no-stress boundary conditions in this subsubsection.In this case, sine functions may be used as eigenfunctions as in standard RB convection, and the z derivative of the base state is ∂ z T B ≈ − cos t.
In (4.31), we take (4.34a,b) to arrive at a quadratic equation for ν λ (t): where To find the stability threshold, we solve the quadratic equation (4.35), set the resulting νλ equal to zero, and find with the maximization over k in the definition of ν in (4.32) leading to k cr = π/ √ 2, the critical wavenumber in steady RB with the same boundary conditions.The value of λ leading to the largest Ra A is λ ≈ 1.319, leading to Ra A ≈ 2891.38, which is exactly what we find numerically, as seen in figure 5.

Heat-flux modulation
Carrying out the linear and nonlinear stability calculations as described in the preceding sections, we are able to find results for linear stability, asymptotic stability and strong global stability for specified ω, Pr and boundary conditions.Results for no-slip conditions on the top and bottom are shown in figures 2(a,c), while no-stress conditions are shown in figures 2(b,d).As usual in RB convection, no-slip conditions lead to a higher critical Rayleigh number.
For linear stability, we find that the critical Rayleigh number always arises from either Im(μ) = 0 or Im(μ) = 1/2, representing synchronous and subharmonic disturbances, respectively.There is a stark contrast between low and high frequencies.At low frequencies, Ra L generally decreases with ω, but in an oscillatory manner as the critical instability switches between various modes of synchronous and subharmonic disturbances, as shown by the discontinuities in the most dangerous wavenumber curves (figures 2c,d).At high frequencies, the critical instability is always subharmonic, and an asymptotic balance is reached with Ra L ω −2 and k L ω −1/2 approaching non-zero constants, as shown in figure 3.For linear stability, no-slip conditions with Pr = 1 give Ra L ω −2 → 22.58, while no-stress conditions give Ra L ω −2 → 12.44.The nonlinear stability thresholds do not appear to reach the same asymptotic relationship between the critical Rayleigh number and the modulation frequency.The nonlinear stability results do not change as radically with frequency as the linear stability results overall, though the threshold for both asymptotic stability and strong global stability does go down at low frequencies.Note that ω ≈ 10 7 for Lake Superior, and Ra ≈ 10 20 at 3 • C (Christopher 2021).This value of ω is above numerically attainable values, which makes the large-ω limit results interesting.Both of these values use molecular thermal diffusivity.As the temperature approaches the temperature of maximum density for water, the coefficient of thermal expansion approaches zero.The Rayleigh number is proportional to the coefficient of thermal expansion, so at some point the Rayleigh number must pass through the critical Rayleigh numbers found here.
The dependence of the critical Rayleigh numbers on Pr is shown in figure 4 for ω = 100.It can be seen that Ra L changes by orders of magnitude as Pr is varied, while Ra A stays in a relatively narrow range.In contrast, strong global stability Ra S is independent of the Prandtl number.A subcritical instability is an instability arising for a Ra value between the linear and nonlinear stability thresholds.There is therefore a very large region for potential subcritical instabilities at low Pr, with the region increasing as Pr decreases, as seen in figure 4. ω ω Low Prandtl number means decreasing viscosity, so on the one hand, it should be easier to trigger instability.However, for small Prandtl number, viscosity on its own disappears from the linear system, which contains Ra Pr.It is possible that both effects compensate for nonlinear stability, leading to a nearly flat curve with a maximum near Pr = 1.It is somewhat surprising, then, to find that Ra A reaches a maximum near Pr = 0.6 and then decreases with decreasing Pr for ω = 100 and no-stress boundary conditions, with a similar result for no-slip conditions.Calculations for ω = 10 show the same pattern of behaviour.
The asymptotic behaviour of Ra L with Pr is readily predicted by looking at the linear equation system (3.6)-(3.7).In the large Prandtl number limit, the first term on the  left-hand side of (3.6) is negligible, so Pr disappears from the linear system, as for the classical RB configuration.Then Ra L is independent of Pr.On the contrary, in the small Prandtl number limit, the second term on the left-hand side of (3.6) is negligible: the only remaining parameter is then Ra L Pr, and the stability threshold in term of this parameter becomes Ra L Pr = const., yielding the scaling Ra L = const.× Pr −1 observed in figure 4(a).

Temperature modulation
Linear stability results for non-zero-mean temperature modulation of one boundary can be found in Or & Kelly (1999).For completeness, we include here nonlinear stability results for that set-up, with Ra defined appropriately for the configuration.These results are shown in figure 5.The general dependence of the critical Rayleigh number is the same as in the modulated flux case, and only the specific numbers are different.
The no-stress case can be treated with sine eigenfunctions, meaning that the stability problem reduces to a single ODE.It is also possible to use the WKB approximation for this case, as in Or ( 2001), but we do not pursue further WKB calculations here for the reasons discussed in § 3.2.
Results scaled for large ω are shown in figure 6.As ω → ∞, the appropriately defined Rayleigh number grows with ω 3/2 , and the critical wavenumber grows with ω 1/2 .For linear stability, no-slip conditions with Pr = 1 lead to Ra L ω −3/2 → 27.86, while no-stress Ra ω -3/2 conditions lead to Ra L ω −3/2 → 18.38.The nonlinear stability thresholds do not appear to reach the same asymptotic relationship between the critical Rayleigh number and the modulation frequency.

High frequency
In this subsection, we look at high-frequency results for all configurations to compare the behaviour of the linear and global stability thresholds.Though we have treated explicitly only the zero-temperature top boundary condition listed in (2.11), it is of course possible to use a no-flux top boundary condition.As the modulation frequency is increased, the base state temperature profile becomes largely confined to a small layer near the modulated surface at the bottom.For large ω, the base state for heat-flux modulation in (2.14) leads to the following form for the base state derivative: (5.1) The (dimensionless) boundary layer thickness is therefore δ = O(ω −1/2 ), so that for large enough ω, we might expect the same results even with different boundary conditions imposed at the top at z = 1 because the base state derivative has hardly any influence there.
Figure 7 shows the critical Rayleigh numbers for all possible configurations.The scaled linear Rayleigh number used in figure 7 in agreement with figures 3 and 6.For all frequencies shown here, the critical disturbance is subharmonic.For linear stability, by ω 100, the boundary conditions at the non-modulated surface have ceased to affect the critical Rayleigh number: Ra ∞ converges towards a constant value that is essentially dependent only on the conditions at the modulated surface.This scaling behaviour can be explained following a local approach similar to that of Howard (1966), assuming that all the dynamics takes place in the boundary layer δ and that the depth d of the system is no longer a relevant parameter of its dynamics.Then one can define a local Rayleigh number as Ra × δ 4 ∼ Ra ω −2 for modulated flux, and Ra × δ 3 ∼ Ra ω −3/2 for modulated temperature.Instability starts once the local Rayleigh number -here the scaled linear Rayleigh number (5.2) -reaches a given critical value that depends only on the conditions at the modulated surface.The most dangerous mode has a wavenumber inversely proportional to the only relevant scale of the system, δ ∼ ω −1/2 , in agreement with figures 3(c,d) and 6(c,d).
In contrast, for nonlinear stability, the non-modulated boundary condition does affect the critical Rayleigh number even at high frequencies.The local analysis in the boundary layer is not relevant.Figure 7 shows that results for the four possible top boundary conditions do not converge at high frequency as they do in linear stability.Asymptotic stability results indicate the same pattern, with the top boundary condition influencing results even at higher frequencies.For nonlinear stability, the critical Rayleigh numbers grow at a rate closer to Ra ∼ ω as ω → ∞.Considering the scaling (5.2) for linear stability, this means that the potential region for subcritical instabilities grows rapidly as ω → ∞.

Validation by direct numerical simulations
Our purpose here is to validate the main features of our analytical stability analysis by performing initial value, two-dimensional DNS of the full equations, starting from the purely diffusive base state solution (2.14) plus some infinitesimal perturbations of the temperature field.A complete numerical study of the system -including, for instance, bistability analysis in the range below the linear threshold and above the nonlinear one, The values of the linear critical Rayleigh number Ra L given here come from the analytical study (see figure 2).
the study of the highly non-linear dynamics at large Rayleigh number, or the processing of more realistic boundary layer forcing -is left for future work.
6.1.Numerical method Equations (2.7)-(2.9)with temperature boundary conditions (2.10)-(2.11)are solved using the commercial software COMSOL Multiphysics, based on the finite element method.Note that for numerical efficiency, it is better to start with a zero flux at the bottom: hence our bottom forcing is ∂ z T = sin t at z = 0, shifted by π/2 compared to the theoretical study, with no further consequences.The computational domain is rectangular, with dimensionless depth 1 and dimensionless width Γ 8, chosen to include at 4 wavelengths of the first excited mode.Top/bottom velocity boundary conditions are either no-slip or no-stress, and we use periodic boundary conditions in the horizontal direction for both temperature and velocity.The mesh is triangular in the bulk and rectangular close to the top and bottom plates, where it is strongly refined.We use standard Lagrange elements, quadratic for the pressure, and cubic for the velocity and temperature fields.The total number of degrees of freedom is at least 2 × 10 5 .Grid convergence and influence of the aspect ratio Γ were tested for each studied value of the Rayleigh number and forcing frequency.No stabilization technique is used.The implicit, time-dependent solver employs the backward differentiation formula with accuracy order 2-3 and relative tolerance 5 × 10 −3 .We impose a minimum of 50 time steps per forcing period.A random noise of typical amplitude 10 −6 is added to the diffusive temperature field (2.14) as the initial condition.Then the code is run for at least 1.5 diffusive time, or until saturation of the kinetic energy for the unstable configurations.Table 1 lists the characteristics of the 15 simulations used in the results presented in figures 8 to 12.Many other simulations were performed to confirm the trends shown here, but are not presented.

Linear and nonlinear stability
We first checked the linear stability results.To do so, we performed a number of DNS runs, systematically changing the Rayleigh number around the theoretical critical value Ra L determined in § 5. We then plot the space-averaged value of the kinetic energy as a function of time: after a short transient due to the adjustment of the initial, random perturbation of the temperature field, it is well fitted by an exponential function of type K 0 exp [2σ (t − t 0 )], with σ the instability growth rate (see e.g.figure 8a).An example of a systematic study for ω = 100 and Pr = 1.0 is shown in figure 8(b).The threshold for instability (where σ = 0) is in perfect agreement with the theoretical prediction.
Our numerical code with a random, infinitesimal, initial temperature perturbation is not well fitted to study the nonlinear stability, which would require imposing as the initial condition the most dangerous mode in all the velocity, pressure and temperature fields.We can nevertheless check the existence of the different regimes.Figure 9 shows the space-averaged kinetic energy for ω = 100, Pr = 1.0, and four different Rayleigh numbers: just above the linear threshold, Ra = 1.025RaL , just below it, Ra = 0.95Ra L , just below the asymptotic nonlinear threshold, Ra = 0.95Ra A , and just below the strong nonlinear threshold, Ra = 0.95Ra S .The main expected features of the different regimes are recovered: above the linear threshold, the small perturbation grows exponentially in time, while below the strong nonlinear threshold, it decreases exponentially.In between, the disturbance energy might grow transiently during a cycle, but for the infinitesimal initial perturbations considered here, it always experiences overall net decay.Again, this is not a complete study of the nonlinear stability, which would require more advanced DNS, but it illustrates the sufficient conditions provided by the nonlinear stability results.

Synchronous and subharmonic modes
One of the most surprising results from the linear analysis is the competition between synchronous and subharmonic modes at a relatively low forcing frequency ω.To verify this, we have performed simulations for various ω, just above the stability threshold.Results are shown in figure 10 and confirm the analytical results.Note that the mode selection is very sensitive to the aspect ratio Γ because of the influence of the imposed periodicity on the wavelength selection.For instance, convergence of the results shown in figure 10 was not reached for Γ = 8.0 used in the previous DNS.
Figures 11 and 12 allow us to further understand the origin of the two different modes.The synchronous mode is the most straightforward to understand: heating the system from below leads to a transient destabilization of the otherwise stably stratified system, and instability appears with a period similar to the forcing; negative flux then restabilizes the system, before a new cycle begins.However, this synchronous mode is clearly subdominant close to linear threshold, where most of time a subharmonic mode kicks in first.From figure 12, both modes correspond to a similar velocity pattern, i.e. one big cell over the whole depth.This cell is mostly stationary, but the direction in which the fluid flows along this cell might reverse or not between two successive forcing cycles, respectively leading to subharmonic and synchronous modes.(Note also the positive/negative reflection symmetry of the subharmonic signals in figure 11a.)If we look at the temperature field at the end of the decreasing flux part of the first cycle (t/2π = 4.25), then we can notice for the subharmonic case negative temperature perturbations on the left and right sides of the zoom, as opposed to the synchronous case: this might lead to a locally stronger bottom temperature gradient at these locations, hence triggering a rising convective velocity there, and a sinking return flow at the central location, which was formerly rising.This mechanism triggers the instability with a period twice that of the forcing.This process is all the more efficient for large ω, i.e. when the temperature perturbation from the previous cycle does not have time to diffuse away, hence the predominance of subharmonic modes at large ω.Nevertheless, preliminary studies when increasing Ra show that these subharmonic modes are restricted to the close neighbourhood of the stability curve: as soon as buoyancy forcing becomes strong enough, the boundary layer rapidly becomes unstable at each cycle before the building up of any subharmonic interaction, and the readily expected synchronous regime appears.As an illustration, for the case ω = 100 and Pr = 1.0 studied in figure 8, a synchronous regime dominates at Ra = 3.5Ra L at saturation, while the subharmonic mode still persists at Ra = 2Ra L .We expect that the competition between the fine tuning necessary to trigger a subharmonic mode, and the most direct, but less efficient excitation of a synchronous mode, also explains the mode alternation observed at low ω (see e.g.figure 10).non-monotonically as ω decreases.The linear stability problem becomes ill-conditioned for an O(1) frequency.For large enough ω, the critical instability is always subharmonic, and Ra L ω −2 and Ra L ω −3/2 approach an O(10) constant in the modulated flux and modulated temperature cases, respectively.The critical Rayleigh numbers for nonlinear stability grow more slowly with ω, approximately linearly.The nonlinear stability results complement the linear stability results, showing that the window of possible Rayleigh numbers for subcritical instability is relatively small for low frequencies but increases rapidly as ω → ∞.
The modulated flux set-up considered here is relevant to situations in nature where a body of fluid experiences periodic heating at the surface, such as the diurnal heating of a lake by the sun.The model in this paper uses a zero-mean heat-flux modulation profile at the boundary, meaning that the net heat flux over each period is zero.This is a simplification of the motivating case of springtime warming of ice-free Lake Superior because the lake warms up during the spring.Despite this difference, the results in this paper may provide insight at the time when the Rayleigh number passes from supercritical to subcritical as the coefficient of thermal expansion goes from negative to positive.The most realistic boundary conditions for the lake would be modulated heat flux and no-stress conditions at the free surface, and zero heat flux and no-slip conditions at the lake bottom.Another example of this set-up arising in the analysis of natural phenomena is Coenen et al. (2021).
We have treated the modulated flux condition at one boundary as being representative of radiative heating confined to a thin layer near the surface, and have also neglected effects from rotation.Future work could include these additional factors.Radiatively driven convection without modulation has recently been used experimentally in Bouillaut et al. (2019) to observe the transition to the ultimate scaling regime of RB convection, where the Nusselt number scales with the square root of the Rayleigh number.Radiative heating could be incorporated into the stability methods used here, and the theoretical modulation profile would then need to avoid radiative cooling.
When considering linear stability, rotation generally has a stabilizing effect on RB convection, as shown in Chandrasekhar (1961), and we would expect the same effect when combined with modulation.When considering nonlinear stability, the form of the energy used here in the energy method is not sensitive enough to include rotation because the inner product of the velocity with the Coriolis term is zero.To find nonlinear stability results with rotation, researchers have had to use a modified energy that leads only to conditional stability results, as detailed in Galdi & Straughan (1985), for example.

Figure 3 .Figure 4 .
Figure 3.As figure 2, but with results now scaled for large ω and computed over a larger range of ω.

Figure 6 .
Figure 6.As figure 5, but with results now scaled for large ω and computed over a larger range of ω.

Figure 7 .
Figure 7.Comparison of results with Pr = 1 for all 16 possible combinations of boundary conditions and modulation style: no-slip or no-stress for velocity, zero-temperature or no-flux for temperature, and heat-flux or temperature modulation.(a) Linear stability; (b) global stability.Colour indicates the velocity boundary condition at the surface of modulation and the modulation type, with the associations listed in the legend.Line style indicates the boundary conditions at the non-modulated surface: solid indicates no-slip and zero-temperature; dashed indicates no-slip and no-flux; dotted indicates no-stress and no-flux; dash-dotted indicates no-stress and zero-temperature.

961Figure 8 .
Figure 8. Linear stability study for ω = 100 and Pr = 1.0.(a) Two examples of the time evolution of the space-averaged kinetic energy and of the determination of the exponential growth rate.(b) Measured growth rate as a function of the Rayleigh number for 9 runs.A complete list of parameters is provided in table 1.

961Figure 9 .Figure 10 .
Figure 9. Temporal evolution of the space-averaged kinetic energy as a function of time for four DNS runs illustrating the linear and nonlinear stability regimes at ω = 100, Pr = 1.0, and for Ra = 1.025Ra L , 0.95 Ra L , 0.95 Ra A and 0.95 Ra S , respectively.Values of Ra L , Ra A and Ra S come from the analytical study (figure 3).A complete list of parameters is provided in table 1.

Figure 11
Figure 11.(a) Time evolution over two forcing periods of the imposed bottom heat flux, of the space-averaged kinetic energy, and of the perturbation temperature and vertical velocity at the centre of a 'hot' cell close to the middle of the computational domain, i.e. at x = 8.4, z = 0.2 for the synchronous case (left) and x = 7.6, z = 0.2 for the subharmonic case (right).The three variables are rescaled to appear on the same y-axis.(b) Snapshot at time t/2π = 4 of the perturbation temperature field normalized by the maximum value over the two cycles shown in (a), and of the streamlines of the associated field.The stars show the locations where the local data in (a) are taken.A complete list of DNS parameters is provided in table 1.

Table 1 .
Dimensionless numbers and velocity boundary conditions for the 15 simulations used in figures 8-12.