Ultimate heat transfer in `wall-bounded' convective turbulence

Direct numerical simulations have been performed for turbulent thermal convection between horizontal no-slip, permeable walls with a distance $H$ and a constant temperature difference $\Delta T$ at the Rayleigh number $Ra=3\times10^{3}-10^{10}$. On the no-slip wall surfaces $z=0$, $H$ the wall-normal (vertical) transpiration velocity is assumed to be proportional to the local pressure fluctuation, i.e. $w=-\beta p'/\rho, +\beta p'/\rho$ (Jim\'enez et al., J. Fluid Mech., vol. 442, 2001, pp. 89-117), and the property of the permeable wall is given by the permeability parameter $\beta U$ normalised with the buoyancy-induced terminal velocity $U={(g\alpha\Delta TH)}^{1/2}$, where $\rho$, $g$ and $\alpha$ are mass density, acceleration due to gravity and volumetric thermal expansivity, respectively. A zero net mass flux through the wall is instantaneously ensured, and thermal convection is driven only by buoyancy without any additional energy inputs. The critical transition of heat transfer in convective turbulence has been found between the two $Ra$ regimes for fixed $\beta U=3$ and fixed Prandtl number $Pr=1$. In the subcritical regime at lower $Ra$ the Nusselt number $Nu$ scales with $Ra$ as $Nu\sim Ra^{1/3}$, as commonly observed in turbulent Rayleigh-B\'enard convection. In the supercritical regime at higher $Ra$, on the other hand, the ultimate scaling $Nu\sim Ra^{1/2}$ is achieved, meaning that the wall-to-wall heat flux scales with $U\Delta T$ independent of the thermal diffusivity, although the heat transfer on the wall is dominated by thermal conduction. The physical mechanisms of the achievement of the ultimate heat transfer are presented.


Introduction
The flow driven by buoyancy is called thermal convection, and it plays an important role in a wide variety of phenomena of geophysics, astrophysics and engineering applications. One of the canonical configurations of thermal convection is the Rayleigh-Bénard convection (RBC) observed in a horizontal fluid layer heated from below and cooled from above. In RBC, buoyancy forcing is characterised in terms of the Rayleigh number Ra and the flow becomes turbulent eventually as Ra increases.
It is known that the Nusselt number N u (dimensionless vertical heat flux) exhibits the power law of Ra, N u ∼ Ra γ , for a certain value of γ in the turbulent state of RBC. For more than half a century, various predictions have been made to clarify the scaling exponent γ. Priestley (1954) derived γ = 1/3 from similarity argument, and Malkus (1954) also led to γ = 1/3 based on the assumption that heat transfer is determined by the marginal instability of a thermal boundary layer. Kraichnan (1962) predicted γ = 1/2 with a logarithmic correction, N u ∼ P r 1/2 Ra 1/2 (ln Ra) −3/2 , as a scaling in a high-Ra asymptotic state with turbulent boundary layers. The scaling N u ∼ P r 1/2 Ra 1/2 is currently known as the ultimate scaling. It has been derived as a rigorous upper bound on the heat transfer in RBC by applying variational methods to the Boussinesq equations (Doering & Constantin 1992, 1996Plasting & Kerswell 2003), and has recently been obtained as a maximal heat transfer scaling between two parallel plates (Motoki et al. 2018). The ultimate scaling relates to the Taylor's energy dissipation law of high-Reynolds-number turbulence via the rigorous energy budget equation of thermal convection. In the ultimate heat transfer the energy dissipation and the scalar dissipation (corresponding to the vertical heat flux) are independent of the kinematic viscosity or the thermal diffusivity.
Recently, Grossmann & Lohse (2000, 2002, 2011 have proposed a unifying scaling theory in RBC, and its validity has been demonstrated by a lot of experimental and numerical studies (see Ahlers et al. 2009). Their theory is based on the energy budget equation relating the energy and scalar dissipation rates and on the decomposition of the flow field into a boundary layer and a bulk region. The theory gives different scaling laws depending on whether the total energy and scalar dissipation rates are dominated by the bulk or the boundary layer. At a high-Ra regime, in which the contribution from the bulk is dominant, the scaling N u ∼ Ra 1/3 is given if the thermal boundary layer is thinner than the velocity boundary layer, but the ultimate scaling N u ∼ P r 1/2 Ra 1/2 is anticipated if the thermal boundary layer is thicker than the velocity boundary layer.
The question of whether or not the ultimate scaling can be achieved has long attracted a great deal of attention, and much effort has been spent on both experimental and numerical studies in the past few decades. However, N u ∼ P r 1/2 Ra 1/2 has not been found as yet in conventional RBC, and what has been observed experimentally and numerically at high Ra is the scaling N u ∼ Ra 1/3 (see Ahlers et al. 2009;Chillà & Schumacher 2012).
It is known that the ultimate scaling N u ∼ P r 1/2 Ra 1/2 can be observed in turbulent thermal convection without horizontal bounding walls on which thermal and velocity boundary layers should have appeared. Such wall-less thermal convection was numerically examined in a triply-periodic domain with a constant temperature gradient in the vertical direction (Calzavarini et al. 2005), and was experimentally investigated in a vertical tube connecting high-and low-temperature chambers (Pawar & Arakeri 2016). The ultimate scaling has also been reported for the thermal convection in a cylindrical container radiatively heated from below, instead of conventional RBC heating (Lepot et al. 2018;Bouillaut et al. 2019). In the radiatively-driven convection N u ∼ Ra 1/3 has been observed when the thickness of heating layer is thin, and the scaling has been found to change to N u ∼ P r 1/2 Ra 1/2 with increasing the thickness.
In case of conventional RBC heating, it has been found that surface roughness on horizontal walls transiently yields the ultimate scaling N u ∼ P r 1/2 Ra 1/2 in the limited range of Ra where the thermal conduction layer thickness is comparable to the size of roughness elements (Zhu et al. 2017(Zhu et al. , 2019Tummers & Steunebrink 2019). This transient scaling would not imply the transition to the asymptotic ultimate scaling, because a further increase in Ra leads to saturation down to the usual scaling N u ∼ Ra 1/3 . It is still an open question whether or not the ultimate heat transfer can be achieved by introducing an ingenious contrivance, such as wall roughness and so on, into wall-bounded RBC heated conventionally.
In this study, we introduce wall permeability into RBC. Jiménez et al. (2001) have investigated turbulent momentum transfer in numerically simulated porous channel flow to find out that the wall permeability significantly enhances momentum transfer. In their simulation the fluid crosses the porous wall surface with a wall-normal velocity proportional to pressure fluctuations. This boundary condition mimics the behaviour of a zero-pressure-gradient boundary layer over a Darcy-type porous wall (Batchelor 1967, pp. 223-224) with a constant-pressure plenum chamber underneath. We perform direct numerical simulations (DNS) for convective turbulence between horizontal no-slip, massneutral permeable walls with a constant temperature difference for fixed Prandtl number P r = 1 by using Jiménez et al.'s (2001) boundary condition on a permeable wall. We report that the wall permeability brings about the ultimate heat transfer N u ∼ Ra 1/2 at a high Rayleigh number in spite of the presence of a thermal conduction layer on the walls. We inspect scaling laws and turbulence structure in thermal convection between the permeable walls as well as impermeable walls to discuss why the ultimate heat transfer can be achieved by the introduction of permeable walls. This paper is organised as follows. The numerical procedure to solve the Boussinesq equations with the no-slip, permeable boundary conditions is presented in §2, and it is confirmed that there are no additional energy inputs except for buoyancy power in §3. Scaling properties and turbulence structure in thermal convection between permeable and impermeable walls are presented in §4, and the physical interpretation of the scaling laws is provided in §5. The summary and outlook are given in §6. The Prandtl number dependence of the scaling of N u with Ra is briefly shown in appendix A, where it is demonstrated that the ultimate scaling can also be observed for the Prandtl number P r = 7.

Direct numerical simulation
We conduct DNS for turbulent thermal convection between horizontal plates with a distance H and a constant temperature difference ∆T . The Oberbeck-Boussinesq approximation is employed, wherein density variations are taken into account only in the buoyancy term. The two horizontal and vertical directions are respectively denoted by x, y and z (or x 1 , x 2 and x 3 ). The corresponding components of the velocity u(x, t) are given by u, v and w (or u 1 , u 2 and u 3 ) in the horizontal and vertical directions, respectively.
The governing equations are the Boussinesq equations where p(x, t) is the pressure, T (x, t) is the temperature, and ρ, ν, g, α and κ are mass density, kinematic viscosity, acceleration due to gravity, a volumetric expansion coefficient and thermal diffusivity, respectively. e z is a unit vector in the vertical direction. The velocity and temperature fields are supposed to be periodic in the horizontal (x-and y-) directions, and the periods in the x-and y-directions are taken to be L. We suppose that the two horizontal walls are composed of porous media with constantpressure plenum chambers underneath and overhead. The lower (or upper) wall and the associated plenum chamber are heated from below (or cooled from above). On the permeable wall surface the vertical velocity w is assumed to be proportional to the local pressure fluctuation p ′ (Jiménez et al. 2001). The boundary conditions imposed on the walls are where β ( 0) represents the property of permeability, and the impermeability conditions w(z = 0, H) = 0 are recovered for β = 0, while β → ∞ implies zero pressure fluctuations and an unconstrained vertical velocity. The flow situation observed in the thermal convection without horizontal walls (Calzavarini et al. 2005) is intuitively similar to this limit, although not identical. Note that a zero net mass flux through the permeable wall is instantaneously ensured because the transpiration velocity is proportional to the pressure fluctuation with zero mean. We anticipate the no-slip and permeable conditions (2.4) and (2.5) on a wall with a large number of wall-normal through holes. The proportionality coefficient β has the dimension of an inverse velocity, and thus βU represents a dimensionless parameter determining the property of permeable walls if the buoyancy-induced terminal velocity U = (gα∆T H) 1/2 is a proper velocity scale. If the proper velocity scale (say, U w ) is smaller than U as in the subcritical permeable case discussed later (see (5.2) in §5), then the permeable condition βU = const. (= β ′ U w ) to be employed here implies a more permeable wall of larger β ′ (= (U/U w )β). Thermal convection between permeable walls is characterised in terms of the Rayleigh number Ra, the Prandtl number P r and the permeability βU , where (2.7) The vertical heat flux from the bottom to the top wall is quantified by the Nusselt number N u written as where · xyt represents the average over the two horizontal directions and time, and · xyzt is the volume and time average. The rightmost equality is given by the volume and time average of the energy equation (2.3). Let us note that since the walls are isothermal in the permeable and impermeable cases, the temperature fluctuation and so the convective heat flux T w xyt are null on the walls (z = 0, H) at any cases. In the near-wall region, therefore, the conduction heat transfer dominates the convective one even in the permeable case. The Boussinesq equations (2.1)-(2.3) are discretised by employing a spectral Galerkin method based on the Fourier series expansion in the periodic horizontal directions and the Chebyshev polynomial expansion in the vertical direction. The nonlinear terms are evaluated using a spectral collocation method. Aliasing errors are removed with the aid of the 2/3 rule for the Fourier transform and the 1/2 rule for the Chebyshev transform. Time advancement is performed with the third-order Runge-Kutta scheme (or the second-order Adams-Bashforth scheme) for the nonlinear and buoyancy terms and the implicit Euler scheme (or the Crank-Nicolson scheme) for the diffusion terms in the permeable (or impermeable) case. The numerical procedure developed by Jiménez et al. (2001) is applied to satisfy the permeable boundary conditions.
In this paper, we shall present the results obtained from DNS for thermal convection in the impermeable case βU = 0 at Ra = 10 6 − 10 11 and in the permeable case βU = 3 at Ra = 3 × 10 3 − 10 10 for fixed Prandtl number P r = 1 and for fixed horizontal period L/H = 1. The dependence on the Prandtl number is shown in appendix A. We have examined the dependence of heat transfer on the horizontal period in the range of 1 L/H 4 to confirm that the ultimate scaling N u ∼ Ra 1/2 to be shown in §4 can also be achieved for smaller βU in a wider periodic box of larger L/H.

Energy budget
In this section, we discuss the total energy budget in thermal convection between noslip, permeable walls. By taking the volume and time average of an inner product of the Navier-Stokes equation (2.2) with the velocity u and taking account of the boundary conditions (2.4) and (2.5), we obtain is a total energy dissipation rate per unit mass. The left-hand side of (3.1) represents buoyancy power (energy input), while the second and the third terms on the righthand side denote pressure power on the permeable walls and outflow kinetic energy across the permeable walls, respectively. The pressure power on the permeable walls is strictly greater than zero, so that it is always an energy sink. Although its sign cannot be specified rigorously, we have confirmed numerically that the outflow kinetic energy across the permeable walls is also positive in the present DNS, implying that the kinetic energy flows out of the system across the permeable walls. It turns out that as in the impermeable case, thermal convection between the permeable walls is sustained only by the buoyancy power without any additional energy inputs. It has also been found numerically that the pressure power is comparable with the energy dissipation whereas the outflow kinetic energy is much less than the dissipation. The energy to be lost in the system via the permeable walls could be considered to be supplied to the other system, i.e. the flow in porous media, to eventually dissipate therein. The rightmost equality of (2.8) yields the relation among the buoyancy power, the Prandtl number, the Rayleigh number and the Nusselt number given by Substituting (3.3) into (3.1) and taking into account the flow symmetries, we arrive at . The red and blue lines indicate the ultimate scaling N u ∼ Ra 1/2 and the ordinary scaling N u ∼ Ra 1/3 , respectively. The inset shows N u compensated by Ra 1/2 in the permeable case.
Note that in the impermeable case, i.e. conventional RBC, the energy budget is given by (3.5)

Scaling properties and turbulence structure
Let us first discuss the scaling property of the Nusselt number N u with the Rayleigh number Ra. Figure 1 shows N u as a function of Ra. It can be seen that the wall permeability leads to significant heat transfer enhancement over the entire Ra range. In the impermeable case βU = 0 the present DNS data in the horizontally-periodic domain are good agreement with the turbulent data obtained from the experiments (Chavanne et al. 2001;Niemela & Sreenivasan 2006) and the numerical simulation (Stevens et al. 2010) performed in cylindrical containers. At high Rayleigh number Ra ∼ 10 8 -10 10 , N u can be seen to scale with Ra as N u ∼ Ra 1/3 , nearly consistent with the well-known turbulence scaling (see e.g. He et al. 2012). In the permeable case βU = 3, on the other hand, the ultimate scaling N u ∼ Ra 1/2 can be observed at higher Rayleigh number Ra ∼ 10 7 -10 10 , whereas the ordinary scaling N u ∼ Ra 1/3 is confirmed at lower Rayleigh number Ra ∼ 10 6 -10 7 . It is worthy to note that the scaling property of N u critically changes around Ra ∼ 10 7 from N u ∼ Ra 1/3 to N u ∼ Ra 1/2 with increasing Ra.
This critical transition in the permeable case can also be confirmed undoubtedly for the root-mean-square (RMS) vertical velocity w rms = w 2 1/2 xyt on the wall as shown in figure 2. In the subcritical Ra range 10 6 Ra 10 7 , the wall-normal transpiration velocity is   weak in the sense that it is of the order of Ra −1/6 U (see figure 2a), corresponding to the vertical velocity scale in the near-wall region of RBC for P r ∼ 1, i.e. the impermeable case, in which the ordinary scaling N u ∼ Ra 1/3 has been observed. In the supercritical Ra range 10 7 Ra 10 10 , on the other hand, the RMS velocity on the wall is significantly strong in the sense that it scales with the buoyancy-induced terminal velocity U (see figure 2b). In §5, for the case of P r ∼ 1, the near-wall vertical velocity scale Ra −1/6 U will be related with the ordinary scaling N u ∼ Ra 1/3 , and the relevance of the vertical velocity scale U to the ultimate scaling N u ∼ Ra 1/2 will also be discussed.
We would like to stress that the ultimate heat transfer is not simplistically a consequence of just the wall permeability. As will be shown later in this section, the wall permeability can trigger a critical change in convection states, consequently leading to the ultimate scaling N u ∼ Ra 1/2 .
Next we differentiate mean temperature profiles between the supercritical permeable case βU = 3 at Ra 10 7 and the impermeable case βU = 0. Figure 3 presents the mean temperature profiles in the impermeable and permeable cases. In the impermeable case, at higher Ra the profile becomes flatter in the bulk region, while the near-wall temperature gradient becomes steeper. In short the mean temperature profile T xyt /∆T cannot scale with z/H. The behaviour of the mean temperature in the subcritical case at Ra 10 7 similar to that in the impermeable case. In contrast to the impermeable case and the subcritical case, the mean temperature profile in the bulk region seems to scale with ∆T as a function of z/H in the supercritical permeable case at Ra 10 7 , and there remains a finite value of the temperature gradient, i.e., the order of ∆T /H, therein even at high Ra. This contrast should be a crucial consequence of the ultimate heat transfer as will be discussed in §5.
In the permeable case with isothermal wall boundaries, different from the thermal convection without horizontal walls (Calzavarini et al. 2005;Pawar & Arakeri 2016), there exists a thermal conduction layer on the wall, where heat transfer by conduction dominates over that by convection. In figure 3(c,d ) are shown the mean temperature profiles 1 − T xyt /∆T as a function of z/δ, where δ is the thickness of a thermal conduction layer defined as All the profiles in the impermeable case collapse onto a single curve in the thermal conduction layer z/δ 1. It is also the case in the permeable case; however, the thermal conduction cannot be dominant around z/δ ∼ 1 where the convection is also important. The large difference of the temperature profiles at z/δ 1 in the supercritical permeable case at Ra 10 7 implies its reasonable scaling with z/H, shown in figure 3b. As mentioned before, the vertical velocity fluctuation on the permeable walls scales with Ra −1/6 U at subcritical Rayleigh number Ra 10 7 . In the impermeable case (in addition to the subcritical permeable case) the near-wall RMS vertical velocity w rms = w 2 1/2 xyt also scales with Ra −1/6 U as a function of z/δ (see figure 4a). However, the vertical velocity fluctuation in the supercritical case at Ra 10 7 exhibits quite distinct behaviour from that in the impermeable and subcritical cases (see figure 4b).
Figures 5(a-c,d-f) show the RMS vertical velocity normalised by the velocity scale Ra −1/18 U and the buoyancy-induced terminal velocity U , respectively. In the impermeable case at 10 6 Ra 10 11 and the subcritical permeable case at 10 5.6 Ra 10 6.8 , the RMS vertical velocity in the bulk region is seen to scale with Ra −1/18 U corresponding to the vertical velocity scale in the bulk region of RBC for P r ∼ 1 ( figure 5a,b), and thus it decreases relatively with respect to U as Ra increases (figure 5d,e). In the supercritical permeable case at 10 7 Ra 10 10 , on the other hand, the velocity fluctuation in the bulk is found to scale with U (figure 5f). The near-wall gradient of w rms with respect to z/H in figure 5 is steeper at higher Ra in the impermeable and the subcritical permeable cases, but the same is not true of the supercritical permeable case. Although w rms is not null on the permeable walls as already shown in figure 2, the ratio of near-wall w rms to bulk w rms should be of the order of Ra −1/9 in the subcritical case, implying that the near-wall vertical velocity fluctuation becomes smaller than that in the bulk region at The RMS temperature T rms = (T − T xyt ) 2 1/2 xyt normalised by the temperature scale Ra −1/9 ∆T and the temperature difference ∆T between the walls is shown in figures 6(a-c,d-f), respectively. In the bulk region of the impermeable and subcritical permeable cases, the RMS temperature is seen to scale with Ra −1/9 ∆T ( figure 6a,b), and so it decreases as Ra increases. On the other hand, the temperature fluctuation in the supercritical permeable case is found to scale with ∆T ( figure 6f). This remarkable difference in the scalings of the temperature fluctuation originates from the scaling difference in the mean temperature (cf. figure 3). In the supercritical permeable case the vertical fluid motion across the sustaining mean temperature difference of O(∆T ) in the bulk region can induce the temperature fluctuation of O(∆T ) even at higher Ra, but in the impermeable and the subcritical cases the vanishing mean temperature difference means the small temperature fluctuation. In §5 we shall discuss the different scaling properties of the RMS vertical velocity with Ra −1/18 U and U as well as the difference in scaling of the temperature fluctuation with Ra −1/9 ∆T and ∆T .
Let us now look into turbulence structure of thermal convection. Figure 7 visualises the instantaneous thermal and vortical structures in the impermeable and supercritical permeable case at Ra = 10 9 . The high-temperature thermal plumes are represented by In the impermeable case the small-scale hot plumes are confined to the near-wall region. In contrast, high-temperature plumes of a remarkably large horizontal length scale fully extend from the bottom wall to the top wall through the bulk in the supercritical permeable case, so that heat transfer is highly enhanced. Recently, such promotion of large-scale circulation has been reported for the convective turbulence, which exhibits the ultimate scaling N u ∼ P r 1/2 Ra 1/2 , in the radiatively-driven convection (Lepot et al.

2018) and in the thermal convection between rough walls (Tummers & Steunebrink 2019).
Although the intensity and the size of small-scale tubular vortices playing a role in energy dissipation are different between the impermeable and the supercritical permeable cases, their spatial structure is more or less the same.
In figure 8 are shown the snapshots of the convective heat flux wT (which is proportional to local buoyancy power) on the horizontal plane in the conduction layer height z/δ ≈ 1 and on the midplane z/H = 1/2. Note that in these figures, wT is normalised so that its mean and standard deviation may be zero and unity, respectively, in each plane of the impermeable and the supercritical permeable cases. The spatial distribution near the wall differs greatly between the impermeable and supercritical permeable cases ( figure 8a,b). The near-wall small-scale structures, corresponding to thermal plumes, can be observed in the impermeable case, while the large-scale structure, which is the part of the fully extended large-scale plume, appears even in the vicinity of the wall in the supercritical case. On the midplane there is no significant difference between the impermeable and supercritical cases ( figure 8c,d ). In the bulk region the heat transfer the supercritical permeable case βU = 3. The heat flux wT on the horizontal plane is normalised so that its mean and standard deviation may be zero and unity, respectively.
is dominated by large-scale convection, regardless of the difference in the near-wall dominant thermal structures.
In figure 9 we show the one-dimensional premultiplied buoyancy-power spectra k y kx P (k x , k y , z) as a function of the distance to the wall, z, and the wavelength in the horizontal (y-) direction, λ = 2π/k y . The buoyancy-power spectra P (k x , k y , z) is given by where (·) represents the Fourier coefficients, (k x , k y ) are the wavenumbers in the horizontal (x-and y-) directions, † denotes the complex conjugate and · t is the time average. The lateral and longitudinal axes of the figures are normalised by the conduction layer thickness δ. P denotes the spectrum of the energy input by buoyancy, and it is also relevant to the spectrum of the convective heat flux shown in figure 8. In the impermeable case we can see significant buoyancy power at small scales in the vicinity of the wall, z/δ ∼ 10 0 , leading to the near-wall thermal plumes (figure 8a), in addition to greater buoyancy power corresponding to the large-scale convection in the bulk region. The near-wall heat flux determined by the marginal instability of the thermal conduction layer gives us the scaling N u ∼ Ra 1/3 widely observed in turbulent RBC (Malkus 1954). The dashed line in figure 9(a) stands for λ = 10z. The spectral ridge is on this line, implying that the energy-inputted horizontal scale is proportional to the distance to the wall. This observation suggests that the convective heat flux exhibits hierarchical selfsimilar structure near the wall. In the subcritical permeable case at Ra = 10 6 (figure 9b), the Rayleigh number is too low for the small-scale plumes of λ ≪ L(= H) to appear in the near-wall region. In the supercritical case (figure 9c) the spectral peak is located at the large horizontal scale λ/δ ∼ 10 3 (λ/L ∼ 1) in the near-wall region z/δ ∼ 10 0 roughly consistent with the wall-normal position of the spectral peak of the small-scale thermal plumes, suggesting that the large-horizontal-scale plume is generated in the near-wall region by buoyancy to fully extend from there to the other wall as observed in figure 7(b). This near-wall large-scale energy input corresponds to the large-scale convective heat flux shown in figure 8(b). As will be discussed in the next section 5, the ultimate heat transfer N u ∼ Ra 1/2 can be attributed to the generation of this long-wavelength (and so intense) thermal mode near the wall. In the bulk region apart from the walls the energy is inputted at the large horizontal length scale in all the impermeable and permeable cases. We note, however, that just in the supercritical case the energy to be inputted at the large horizontal scale in the bulk is smaller than that in the near-wall region.

Physical interpretation of scaling laws
Here we shall discuss the physical mechanisms of the ordinary scaling N u ∼ Ra 1/3 and the ultimate scaling N u ∼ P r 1/2 Ra 1/2 in turbulent thermal convection between impermeable and permeable walls. In the present study the Prandtl number has been set to unity, i.e., P r = 1, and thus in this section we assume that P r ∼ 1 (or ν ∼ κ).
Let us start with the thermal convection in the impermeable case and the subcritical permeable case, where the temperature profile is flatter in the bulk region at higher Ra and thus temperature variation is confined to the near-wall layer of the thickness of O(δ) (see figure 3a,b). The vertical velocity is strictly zero on the impermeable walls.
In the subcritical permeable case, as shown in figure 2 and figure 4 transpiration has not been activated in the near-wall region although the walls are permeable. In both the impermeable and subcritical cases, therefore, the near-wall vertical velocity is small in comparison to that in the bulk region. We now suppose that in the near-wall layer with the thickness δ ′ of O(δ) and the temperature difference of O(∆T ) where the vertical velocity scale U w is small, the effect of viscosity is significant. In the vertical component of the Navier-Stokes equation (2.2), the viscous term is comparable with the advection term and the buoyancy term, that is, in the near-wall region. The balance (5.1) between the viscous, the advection and the buoyancy terms in the equation of motion determines the near-wall velocity and the length scales as 3) (recall that U = (gα∆T H) 1/2 and H are the buoyancy-induced terminal velocity and the wall distance, respectively). In the present DNS we have confirmed that the vertical velocity near the impermeable and subcritical permeable walls scales with Ra −1/6 U (see figure 2a and figure 4). Since the definition (4.1) of the thermal conduction layer thickness implies that δ ′ ∼ δ = (H/2N u), we arrive at the scaling law which has been observed in RBC (i.e., the impermeable case) as well as in the subcritical permeable case (see figure1). The scaling law N u ∼ Ra 1/3 has already been given by the several arguments on similarity (Priestley 1954), the marginal instability (Malkus 1954) and the bulk contribution to energy and scalar dissipation (Grossmann & Lohse 2000).
In the bulk region of the impermeable and the subcritical cases, where the effects of viscosity or thermal conduction are no longer significant, the characteristic length scale is H instead of δ (and δ ′ ), and the temperature difference with respect to the height difference of O(H) and the vertical velocity scale are supposed to be ∆T ′ and U b , respectively. In this region the advection and the buoyancy terms balance each other out in the Navier-Stokes equation as Rewriting the Nusselt number (2.8) as and taking into consideration the dominance of convective heat transfer and the scaling (5.4), we have Equations (5.5) and (5.7) yield the temperature difference and the velocity scale as ∆T ′ ∼ Ra −1/9 ∆T, (5.8) Equation (5.9) means that the Reynolds number for thermal convection is of the order of Ra 4/9 P r −2/3 , being consistent with the Grossmann & Lohse's (2000) scaling based on the energy and scalar dissipation in the bulk region. It has been confirmed that the vertical velocity and the temperature fluctiation scale with Ra −1/18 U and Ra −1/9 ∆T , respectively, in the bulk region of the impermeable and subcritical permeable cases (see figure 5a,b and figure 6a,b). Next we consider the thermal convection between the supercritical permeable walls. In this case intense vertical transpiration is induced even in the vicinity of the wall in contrast to the impermeable and the subcritical cases (see figure 2). Although the thermal conduction layer still exists on the wall, there is no near-wall layer of significant change in the vertical velocity, suggesting that the effect of the viscosity on the vertical velocity is negligible anywhere. The vertical motion should exhibit the length scale comparable with H (see figure 7b), and the corresponding temperature difference is of the order of ∆T even in the bulk region (recall the temperature gradient of O(∆T /H) in figure 3b and the temperature fluctuation of O(∆T ) in figure 6f). Therefore, the balance between the dominant advection and buoyancy terms in the Navier-Stokes equation (2.2) gives us The balance between the buoyancy power, the energy dissipation and the pressure power in the energy budget (3.4), suggests the Taylor's dissipation law (energy dissipation independent of ν) ǫ ∼ U 3 H (5.13) and the ultimate scaling N u ∼ P r 1/2 Ra 1/2 ∼ Ra 1/2 , (5.14) where we have taken account of βU ∼ 1. Alternatively, it follows from (5.6) at z/δ ≫ 1 and (5.11) that The ultimate scaling N u ∼ Ra 1/2 has been suggested by Kraichnan (1962) and by Grossmann & Lohse (2000) as high-Ra asymptotics; however, it has not been observed in conventional RBC as yet. In the present DNS of the supercritical permeable case, the vertical velocity has been seen to scale with U in the whole region (see figure 5f), and it has been confirmed that N u ∼ Ra 1/2 at Ra 10 7 (see figure 1).
In the above discussions we have considered the difference in the vertical length scale of thermal convection, δ and H, in the impermeable (and subcritical permeable) case and the supercritical permeable case and its crucial consequences on the scaling properties of heat transfer. The key to the difference in the vertical length scale is the excitation of transpiration in the near-wall region of the permeable wall. As suggested in figure 9, there should be different convection modes of the instabilities in a thermal conduction layer, one of which is the small-scale thermal plume in the impermeable (and subcritical permeable) case, and the other of which is the large-scale plume extending to the other wall in the supercritical permeable case. The excitation of the near-wall transpiration velocity on the permeable wall could be attributed to the different length of convection instability from that on the impermeable wall. In order to identify the different length of the instability, we have performed the linear stability analysis of a conduction state between the impermeable and permeable walls by conducting DNS in conjunction with the Arnoldi iteration. Although the onset of thermal convection in the conduction state is distinct from that in the thermal conduction layer of turbulent convection, we could expect their qualitative similarity. Figure 10 presents the onset Rayleigh number of thermal convection between impermeable and permeable walls. In the impermeable case we confirm the known lowest value Ra c = 1708 for kH = 3.117 (black symbols). In the permeable case, on the other hand, much larger-scale thermal convection (for much smaller kH) can arise from the instability (colour lines with symbols). If such larger-horizontal-scale thermal plume appears in the thermal conduction layer of convective turbulence, then the plume should also possess a larger vertical length scale to induce the significant vertical velocity. Actually the largehorizontal-scale thermal plume has been observed to extend from the near-wall region to the other wall in turbulent convection on the supercritical permeable walls (see figure 7b and figure 8b). The critical transition to the ultimate heat transfer observed in figure 1 and figure 2 could be a consequence of the exchange of near-wall unstable convection modes on the permeable wall.
In DNS of the permeable case we have increased the horizontal period L in the range of 1 L/H 4 and have observed stronger convection in a wider periodic box of larger L/H for βU = 3. This would be because longer-wavelength convection is more significant as a result of the convection instability (see figure 10). Smaller βU ∼ 10 −1 would, however, lead to an optimal length scale of convection and thus no significant dependence of convective turbulence on the horizontal domain size.

Summary and outlook
We have performed the three-dimensional direct numerical simulation (DNS) of turbulent thermal convection between horizontal no-slip, permeable walls with a distance H and a constant temperature difference ∆T . On the no-slip wall surfaces z = 0, H the vertical transpiration velocity has been assumed to be proportional to the local pressure fluctuation (Jiménez et al. 2001), i.e. w = −βp ′ /ρ, +βp ′ /ρ mimicking a Darcytype permeable wall (Batchelor 1967, pp. 223-224). A zero net mass flux through the permeable wall is instantaneously ensured, and convective turbulence is driven only by buoyancy without any additional energy inputs. The permeability parameter is set to βU = 0 (an impermeable case) and βU = 3 (a permeable case) where U = (gα∆T H) 1/2 is the buoyancy-induced terminal velocity. DNS has been carried out at the Rayleigh number up to Ra = 10 11 in the impermeable case and Ra = 10 10 in the permeable case for fixed Prandtl number P r = 1. We have found that the wall permeability leads to the critical transition of the Nusselt number scaling with the Rayleigh number from N u ∼ Ra 1/3 to the ultimate scaling N u ∼ Ra 1/2 as Ra increases.
In the subcritical regime 10 6 Ra 10 7 we have found the scaling law N u ∼ Ra 1/3 commonly observed in turbulent Rayleigh-Bénard convection (RBC) although on the permeable wall, there are weak vertical velocity fluctuations of the order of Ra −1/6 U comparable with the velocity scale of near-wall small-scale thermal plumes in RBC (i.e. the impermeable case). The mean temperature gradient becomes small in the bulk region as Ra increases, and temperature fluctuations scale with Ra −1/9 ∆T in the bulk.
In the supercritical regime 10 7 Ra 10 10 , on the other hand, the ultimate scaling N u ∼ Ra 1/2 has been found. In this supercritical regime the mean temperature profile exhibits a steeper gradient in the very-near-wall thermal conduction layers at higher Ra while a finite value of the temperature gradient remains in the bulk region, implying temperature fluctuations of O(∆T ), in contrast to the vanishing bulk temperature gradient in the impermeable and subcritical permeable cases. This situation is very different from convective turbulence without horizontal walls (Calzavarini et al. 2005;Pawar & Arakeri 2016), in which there is no thermal conduction layer and the ultimate scaling has also been observed. In the supercritical case the significant transpiration velocity is induced even in the vicinity of the wall. The vertical velocity fluctuation scales with U at any height. Although the vertical velocity fluctuation is suppressed near the permeable wall in comparison to the bulk region, there is no near-wall layer of large change in the vertical velocity, suggesting that the effect of viscosity is negligible even in the near-wall region. In such 'wall-bounded' convective turbulence the vertical fluid motion exhibits the large length scale of O(H) in the whole region, and the buoyancy acceleration by the temperature difference of O(∆T ) can achieve the vertical velocity comparable with the terminal velocity U . The ultimate heat transfer is attributed to the resulting large-scale strong plumes extending from the near-wall region of one permeable wall to the other wall. The balance between buoyancy power, energy dissipation and pressure power on the permeable walls in the total energy budget equation provides us with the Taylor's dissipation law ǫ ∼ U 3 /H as well as the ultimate scaling N u ∼ Ra 1/2 . The key to the achievement of the ultimate heat transfer is the activation of transpiration in the near-wall region of the permeable wall, leading to the large-scale and so intense vertical fluid motion. The excitation of transpiration is considered to be a consequence of near-wall larger-horizontal-scale unstable convection mode on the permeable wall, distinct from that on the impermeable or less permeable wall.
Finally, we would like to suggest the possibility of the ultimate heat transfer in physical experiments. The properties of the present permeable wall can be estimated as a porous wall of many fine through holes in the vertical direction with a constant-pressure plenum chamber underneath (or overhead). We install so many holes in the wall that the entire surface of the wall is almost covered by the holes. Supposing the flow through the holes to be laminar and thus be represented by the Hagen-Poiseuille flow, we have its mean velocity w = d 2 32νl ∆p ρ , (6.1) where d, l and ∆p represent the diameter of the holes, the thickness of the wall and the pressure drop through the wall (or the pressure difference with respect to the constant pressure in the plenum chamber), respectively. From the permeable boundary condition (2.5), the permeability coefficient β can be expressed rigorously as where we have used the reasonable relation w 2 xyt wall ∼ w 2 . It follows from (6.3) that d/H ∼ (βU ) 1/2 P r 1/4 Ra −1/4 . (6.6) Now we map the above estimates onto turbulent thermal convection between the permeable walls. In the supercritical permeable case βU = 3 at Ra ∼ 10 9 for P r ∼ 1, equation (6.6) tells us that d/H ∼ 10 −2 . Since the RMS vertical velocity on the supercritical permeable wall is approximately 10% of U (see figure 2b), the Reynolds number of the flow in the holes might be Re = wd/ν ∼ 10 1 , implying that the flow is laminar. Therefore, we may say that the properties of the supercritical permeable walls can be implemented by using the above porous walls to achieve the ultimate heat transfer in physical experiments.

Appendix A. Prandtl-number dependence
In order to examine the effects of the Prandtl number P r on the scaling of the Nusselt number N u with the Rayleigh number, we have performed DNS of turbulent convection between the permeable walls for P r = 7. We inspect the ultimate scaling law N u ∼ P r 1/2 Ra 1/2 .
(A 1) The filled and open circles respectively represent the present DNS data for P r = 7 and P r = 1 in the permeable case βU = 3. The red and blue line indicate N u ∼ P r 1/2 Ra 1/2 and N u ∼ Ra 1/3 , respectively. The inset shows N u compensated by P r 1/2 Ra 1/2 .  Figure 12: RMS vertical velocity on the wall z = 0 normalised by (a) Ra −1/6 P r 1/6 U and (b) U in the permeable case βU = 3 for P r = 7. Figure 11 shows N u compensated by P r 1/2 as a function of the Rayleigh number Ra at P r = 7 and P r = 1 for the horizontal period L/H = 1 and the permeability βU = 3. The transition to the ultimate scaling N u ∼ P r 1/2 Ra 1/2 is also observed for P r = 7. In the supercritical Ra range, the compensated N u-plots roughly collapse on a single line. The scaling behaviour in the subcritical Ra-range for P r = 7 is different from that for P r = 1, and the transition point seems to have a slight P r dependence. The critical transition is also observed in the RMS vertical velocity on the wall for P r = 7 as shown in figure 12. We have confirmed that the other turbulent statistics and structures are similar to those observed for P r = 1.