2.1. Governing equations

Rayleigh–Bénard (RB) convection is described by the incompressible Navier–Stokes equations coupled to buoyancy effects under the Boussinesq approximation. In non-dimensional form, the equations read

(2.1)$$\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} = \sqrt{\frac{Pr}{Ra}}\nabla^2 \boldsymbol{u} - \boldsymbol{\nabla} p + T\boldsymbol{e}_y, \end{gather}$$

(2.2)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$

(2.3)$$\begin{gather}\frac{\partial T}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} T = \frac{1}{\sqrt{Pr Ra}}\nabla^2 T. \end{gather}$$

We restrict to two spatial dimensions in this work. Here, $\boldsymbol {u}$ denotes velocity, $p$ the pressure and $T$ the temperature. We denote by $u$ and $v$ the horizontal and vertical velocity components, respectively. Buoyancy effects are included in the momentum equation by means of the term $T\boldsymbol {e}_y$, where $\boldsymbol {e}_y$ denotes the vertical unit vector. The dimensionless parameters that determine the flow are the Rayleigh number $Ra = g\beta \Delta L_y^3/(\nu \kappa )$ and the Prandtl number $Pr=\nu /\kappa$. The Rayleigh number describes the ratio between buoyancy effects and viscous effects and is set to $10^{10}$ to set the focus on the challenging high-$Ra$ convection regime. The Prandtl number determines the ratio of characteristic length scales of the velocity and the temperature, and is set to 1. The gravitational acceleration is denoted by $g$, the thermal expansion coefficient by $\beta$, the temperature difference between the boundaries of the domain by $\varDelta$, the kinematic viscosity by $\nu$ and the thermal diffusivity by $\kappa$. The computational domain is rectangular with horizontal length $L_x$ and vertical length $L_y$, which are chosen as 2 and 1, respectively. The domain is periodic for all variables in the horizontal direction and wall-bounded in the vertical direction. No-slip boundary conditions are imposed for the velocity at the walls. The non-dimensional temperature is prescribed at 1 at the bottom wall and 0 at the top wall.

Two-dimensional RB convection is fundamentally different from three-dimensional RB convection. The main advantage of restricting the flow to two spatial dimensions is that this enables DNS at a very large Rayleigh number (Zhu et al. Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018). Additionally, the large-scale circulation that appears in three-dimensional RB convection displays quasi-two-dimensional behaviour and shows strong similarities with the large-scale circulation in two-dimensional RB convection (van der Poel, Stevens & Lohse Reference van der Poel, Stevens and Lohse2013).

2.2. Numerical methods

The adopted spatial discretization is an energy-conserving finite difference method (Vreman Reference Vreman2014) and is parallelized as by Cifani, Kuerten & Geurts (Reference Cifani, Kuerten and Geurts2018). A staggered grid arrangement is used for the velocity, the pressure is defined at the cell centres and the temperature is defined on the same grid as the vertical velocity. The latter choice ensures that no interpolation errors occur when computing the buoyancy term in (2.1) (van der Poel et al. Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015). A uniform grid spacing is used along the horizontal direction whereas a hyperbolic tangent grid profile is adopted along the vertical direction. The non-uniform grid realizes refinement near the walls to resolve the boundary layer.

Time integration is performed using the fractional-step third-order Runge–Kutta (RK3) scheme for explicit terms combined with the Crank–Nicholson (CN) scheme for implicit terms, as presented by van der Poel et al. (Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015). Every time step, from $t^n$ to $t^{n+1}$, consists of three sub-stages, of which the steps are outlined below. The superscript $j$, $j=0, 1, 2,$ denotes the sub-stage, where $j=0$ coincides with the situation at $t^n$. In each stage, a provisional velocity $\boldsymbol {u}^*$ is computed according to

(2.4)\begin{equation} \frac{\boldsymbol{u}^* - \boldsymbol{u}^j}{\Delta t} = \left[ \gamma_j H^j + \rho_j H^{j-1} - \alpha_j\mathcal{G}p^j + \alpha_j\mathcal{A}_y^j \frac{\left(\boldsymbol{u}^* + \boldsymbol{u}^j \right)}{2} \right]. \end{equation}

The parameters $\gamma, \rho$ and $\alpha$ are the Runge–Kutta coefficients, given by $\gamma =[8/15, 5/12, 3/4], \rho =[0, -17/60, -5/12]$ and $\alpha =\gamma +\rho$ (Rai & Moin Reference Rai and Moin1991; van der Poel et al. Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015; Cifani et al. Reference Cifani, Kuerten and Geurts2018). Moreover, $H^j$ comprises the discrete convective terms, the discrete horizontal diffusion terms and the source terms. Here, the only source term is the buoyancy term appearing in the evolution of the vertical velocity. The discrete gradient operator is denoted by $\mathcal {G}$. The discrete vertical diffusion, given by $\mathcal {A}_y$, is treated implicitly. The implicit treatment eliminates the severe viscous stability restriction caused by the use of a non-uniform mesh near the boundary (Kim & Moin Reference Kim and Moin1985). Subsequently, the Poisson equation (2.5) is solved for the pressure to impose the continuity constraint (2.2),

(2.5)\begin{equation} \nabla^2\phi = \frac{1}{\alpha_j\Delta t}\left(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^*\right). \end{equation}

Discretely, this equation takes the form

(2.6)\begin{equation} \mathcal{L}\phi = \frac{1}{\alpha_j\Delta t}\left(\mathcal{D}\boldsymbol{u}^*\right). \end{equation}

Here, $\mathcal {L}$ is the discrete Laplace operator $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\nabla }$ and $\mathcal {D}$ is the discrete divergence operator $\boldsymbol {\nabla }\boldsymbol {\cdot }$. The velocity and pressure are subsequently updated according to

(2.7)\begin{equation} \boldsymbol{u}^{j+1} = \boldsymbol{u}^* - \alpha_j\Delta t\left(\mathcal{G}\phi\right) \end{equation}

and

(2.8)\begin{equation} p^{j+1} = p^j + \phi - \frac{\alpha_j\Delta t}{2 Re}\left(\mathcal{L}\phi \right), \end{equation}

after which the velocity $\boldsymbol {u}^{j+1}$ is divergence-free. Time integration of the energy equation (2.3) follows similarly. The newly obtained velocity is used to generate the energy field $T^{j+1}$ analogously to (2.4).

The convective terms are discretized using the QUICK interpolation scheme (Leonard Reference Leonard1979). The diffusive terms are discretized using a standard second-order finite difference method, for both spatial directions. Similarly, the discrete gradient $\mathcal {G}$, the discrete divergence $\mathcal {D}$ and the discrete Laplacian $\mathcal {L}$ are defined using finite differences.

2.3. Altered dynamics under coarsening

The DNS is carried out on a $4096\times 2048$ grid, which has been shown to be a sufficiently high resolution for the chosen Rayleigh number (Zhu et al. Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018). The coarse-grid numerical simulations will be performed on a $64\times 32$ grid. This coarsening introduces significant discretization errors and does not allow for an accurate resolution of the smaller coherent structures present in the flow. The truncation error of the numerical method on this coarse grid introduces artificial dissipation and additional high-pass smoothing native to the discretization method (Geurts & van der Bos Reference Geurts and van der Bos2005). Figure 1 shows a snapshot of the DNS and the coarse-grid simulation, both in statistically steady states, from which the huge implications of such significant coarsening become apparent.

Figure 1. Snapshots of the (a,c,e,g) DNS ($4096\times 2048$ grid) and (b,d, f,h) coarse no-model simulation ($64\times 32$ grid) in statistically steady states. From panels (a,b) to (g,h), we show temperature, pressure, horizontal velocity and vertical velocity.

The temperature at the mentioned resolutions is shown in figure 1(a,b). The coarsened temperature displays only qualitative large-scale agreement with the DNS temperature, both yielding similar plumes of temperature and similar large-scale circulation. Obviously, the persistent small-scale coherent structures visible in the DNS snapshot are lost on the coarse grid. This loss may also be observed for the pressure, depicted in figure 1(c,d). The high-resolution and low-resolution fields exhibit clear qualitative differences, especially considering the absence of distinct low-pressure regions in the coarse-grid result. A better qualitative agreement between the results at the different resolutions is observed for the velocity components, shown in figure 1(e–h). At low resolution, qualitatively the same flow patterns can be observed as appear in the high-resolution results, albeit with decreased magnitude.

The discrepancies between the velocity and temperature fields at the different resolutions clearly pose the challenge we address in this paper. In the next section, we therefore specify our new forcing approach which largely aims to rectify the observed differences. The extent to which this new approach is successful will be scrutinized in § 4.

3.1. Reconstructing the energy spectra

The momentum equation (2.1) and temperature equation (2.3) can be written as a system of complex ordinary differential equations (ODEs) for the mode coefficients through projection onto a Fourier basis. In what follows, we only describe a spectrum-preserving forcing for the momentum. The same derivation is used to define a forcing for the temperature. We note that a spectral decomposition of the velocity and temperature fields is not straightforward due to the presence of boundaries along one spatial dimension. Therefore, we restrict ourselves to one-dimensional periodic horizontal profiles to ensure that the Fourier series is well defined. This choice enables the model to explicitly identify wall-induced features of the flow. Alternatively, after taking into account the non-uniform grid spacing in the wall-normal direction the velocity field can be decomposed by applying a sine transform in the wall-normal direction or by periodically extending the domain and applying a Fourier transform, but these approaches are not pursued in this study.

For a fixed vertical coordinate $y_l$, the profile $\boldsymbol {u}(x, y_l, t)$ is decomposed into Fourier modes and corresponding complex coefficients $c_{k, l}$, where $k$ denotes the wavenumber. The coefficients satisfy the system of ODEs

(3.1)\begin{equation} \frac{\text{d}c_{k, l}}{\text{d}t} = L\left(\boldsymbol{c}, k, l\right), \quad k=0,\ldots,N_x/2, \end{equation}

where $L$ involves the spectral representation of $\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }$, $\boldsymbol {\nabla }$, $\nabla ^2$ and the source term in (2.1). We have already assumed a finite truncation of the number of Fourier modes in this formulation. The vector $\boldsymbol {c}$ comprises all Fourier coefficients up to the largest resolvable wavenumber.

To arrive at a spectrum-preserving forcing, it is sufficient to consider only the magnitude of the Fourier coefficients $|c_{k, l}|$. These evolve according to a system of ODEs

(3.2)\begin{equation} \frac{\text{d}|c_{k, l}|}{\text{d}t} = L_r(\boldsymbol{c}, k, l),\quad k=0, \ldots, N_x/2, \end{equation}

where $L_r$ is introduced to simplify notation. In the model implementation, the definition and explicit computation of $L$ and $L_r$ is not required. Regarding $|c_{k, l}|$ as a stationary stochastic process, we observe that $\mathbb {E}(|c_{k, l}|^2)$ describes the mean energy content of the $k$th Fourier mode at height $y_l$. By the definition of variance, we find that

(3.3) \begin{equation} \mathbb{E}(|c_{k, l}|^2) = \mathrm{var}(|c_{k, l}|) + \mathbb{E}(|c_{k, l}|)^2. \end{equation}

Thus, it is sufficient to obtain an accurate mean value and variance of $|c_{k, l}|$ to achieve the desired energy content in the $k$th Fourier mode. We aim to reproduce the mean value and variance of $|c_{k, l}|$ by augmenting (3.2) with an Ornstein–Uhlenbeck (OU) process,

(3.4)\begin{equation} \text{d}|c_{k, l}| = L_r(\boldsymbol{c}, k, l)\,\text{d}t + \frac{1}{\tau_{k, l}}(\mu_{k, l} - |c_{k, l}|)\,\text{d}t + \sigma_{k, l}\,\text{d}W_{k,l}^t , \end{equation}

where $\mu _{k, l}$, $\tau _{k, l}$ and $\sigma _{k, l}$ are statistical parameters inferred from a sequence of snapshots of the reference solution. The desired mean value of $|c_{k, l}|$ is given by $\mu _{k, l}$, the term $\sigma _{k,l}\,\text {d}W_{k,l}^t$ serves to match the measured variance. Here, $\text {d}W_{k,l}^t$ is a general stochastic process which can be tailored to the data (Ephrati et al. Reference Ephrati, Cifani, Luesink and Geurts2023b), although the common choice is to let it be a Gaussian process with a variance depending on the time step size (Higham Reference Higham2001). We adopt the latter and include the variance scaling in the definition of $\sigma _{k,l}$. The forcing strength is determined by a time scale $\tau _{k, l}$ and can be specified per $k, l$ separately. Detailed specifications of the adopted values of $\mu _{k,l}, \sigma _{k,l}$ and $\tau _{k,l}$ follow shortly. The combination of the original dynamics and the feedback control, including the stochastic term, arises as the continuous-time limit of the 3DVAR data assimilation algorithm (Blömker et al. Reference Blömker, Law, Stuart and Zygalakis2013). The model assumes that the unresolved dynamics can be accurately represented by independent stochastic processes. Interactions in the vertical direction are included via the fully resolved simulations which are base to the forcing. We will demonstrate that this model is capable of producing accurate results.

The model will be applied as a prediction-correction scheme. In the case of the RK3 scheme adopted in this work, the following steps are performed for each RK sub-stage. First, the provisional velocity $\boldsymbol {u}^*$ is computed as in (2.4). The Fourier coefficients $\tilde {c}_{k, l}$ are computed for this provisional velocity, after which the model is applied in the form of a correction. In terms of the Fourier coefficients, the algorithm reads

(3.5)$$\begin{gather} \tilde{c}_{k, l}^{j+1} = c_{k, l}^j + \Delta t \alpha_j L(\boldsymbol{c}^j, k, l)\quad \text{(time integration)}, \end{gather}$$

(3.6)$$\begin{gather}|c_{k, l}^{j+1}| = |\tilde{c}_{k, l}^{j+1}| + \frac{\alpha_j \Delta t}{\tau_{k, l}}\left(\mu_{k, l} - |\tilde{c}^{j+1}_{k, l}|\right) + \sigma_{k, l} \Delta W_{k, l} \quad \text{(correction)}. \end{gather}$$

Here, (3.5) is a simplified notation of the time-marching method presented in § 2.2 and the superscript $j, j=0,1,2,$ denotes the RK sub-stage. The values of $\Delta W_{k, l}$ are samples from a standard normal distribution, independently drawn for each $k$ and $l$. These are determined for each time step and kept constant throughout the sub-stages comprising the time step. The correction only affects the magnitudes of the Fourier coefficients, the phases of $c^{j+1}_{k, l}$ are the same as those of $\tilde {c}^{j+1}_{k, l}$. Velocity fields are subsequently obtained by applying the inverse Fourier transform to the corrected coefficients $c^{j+1}_{k, l}$. After this procedure, the Poisson equation (2.6) is solved using the newly obtained velocity fields and the remaining steps of the sub-stage are completed. Applying the model before solving the Poisson equation ensures that the flow is incompressible at the end of each RK sub-stage. The entire algorithm, with the exception of solving the Poisson equation, may also be applied to the temperature equation. This prediction-correction algorithm has the additional benefit that the model can be easily implemented into already existing computational methods. Applying the correction once every few time steps instead of every time step is common in data assimilation, where real-time data become available sequentially and may not be available at each time step of the numerical simulation. In the current study, all data are collected before performing a coarse numerical simulation and serve to define a forcing term aiming to counteract coarsening effects. The subsequent application of the correction does not noticeably add to the computational effort and is hence applied at each RK sub-stage.

The correction (3.6) will be referred to as nudging. We distinguish between stochastic nudging, which is described by (3.6), and deterministic nudging, for which the stochastic term in (3.6) is omitted. We define a mean $\mu _{k,l, {stoch}}$ and $\mu _{k,l, {det}}$ for these methods, respectively, and specify these below.

The mean $\mu _{k,l}$ is specified such that the desired energy content is reproduced for small values of $\tau _{k,l}$. The magnitude of the coefficients is fully determined by the model in the limit of small $\tau _{k,l}$ and, as a result, (3.3) can be used to derive the mean $\mu _{k,l}$ in this limit. To attain the desired energy contents when using stochastic nudging, we require that $\mu _{k,l, {stoch}}=\mathbb {E}(|c_{k,l}|)$. In the case of deterministic nudging with small $\tau _{k,l}$, the magnitudes of the coefficients remain constant at the value of $\mu _{k,l, {det}}$. Thus, the variance of $|c_{k, l}|$ is set to zero and we require that $\mu _{k,l,{det}}=\sqrt {\mathbb {E}(|c_{k, l}|^2)}$. We note that the limit of small $\tau _{k,l}$ is used to derive the model parameters, but the actual measured values of $\tau _{k,l}$ from the high-resolution data will be (considerably) larger. As a result, the magnitudes obtained using the model will be a combination of the coarse discretization and the high-resolution data.

Treating the evolution operator $L_r$ of the magnitudes of the Fourier coefficients in (3.2) as the identity operator allows us to specify the noise magnitude $\sigma _{k,l}$. This assumption is used in the 3DVAR data assimilation algorithm (Lorenc & Rawlins Reference Lorenc and Rawlins2005) and is valid in statistically stationary states and on short assimilation intervals. This allows replacing $|\tilde {c}^{j+1}_{k,l}|$ by $|c^j_{k,l}|$ and reduces (3.6) to a first-order autoregressive (AR$(1)$) process. We observe that the drift coefficient is $(1-\alpha _j\Delta t/\tau _{k,l})$ and assume that the sample variance $s^2_{k, l}$ is known from the high-fidelity data for every $k, l$. The noise magnitude follows by matching the variance of the AR(1) process with the sample variance, leading to the expression

(3.7)\begin{equation} \sigma_{k,l} = s_{k,l}\sqrt{1 - \left(1 - \frac{\alpha_j\Delta t}{\tau_{k,l}}\right)^2}. \end{equation}

In this paper, the time scale $\tau _{k,l}$ will be defined as the correlation time of the corresponding Fourier coefficient, measured from the high-resolution data. We note that, with the adopted definitions of $\mu _{k,l}$ and $\sigma _{k,l}$, the time scale $\tau _{k, l}$ can take on a range of values whilst still yielding accurate energy spectra. Robustness of the model under variations of $\tau _{k, l}$ is studied in § 4.5. In fact, $\tau _{k, l}$ can take on any positive value larger than or equal to $\alpha _j \Delta t$. Small length scales are expected to yield a small value of $\tau _{k,l}$, resulting in an increased weight towards the model term and an increased noise magnitude for the stochastic term in the nudging. The model term will have a decreased weight at scales for which a large $\tau _{k,l}$ is measured, which is often observed for large spatial scales. These would correspondingly follow the deterministic resolved dynamics more closely.

The proposed prediction-correction method is of the same form as Fourier domain Kalman filtering (Harlim & Majda Reference Harlim and Majda2008; Majda & Harlim Reference Majda and Harlim2012). By defining a prediction and an observation, the approach can be understood as a steady-state filter with a prescribed gain factor and can be placed in the context of data assimilation. The only necessary parameters are the means $\mu _{k,l}$, the variances $\sigma ^2_{k,l}$ and the correlation times $\tau _{k,l}$, which are interpreted as follows. At each sub-stage of the RK3 scheme, the prediction is obtained by evolving the velocity and temperature fields according to the coarse-grid discretization. The ‘observation’ then consists of velocity or temperature fields sampled from the reference statistically stationary state. For stochastic nudging, these are velocity or temperature fields where the magnitudes of the Fourier coefficients are drawn from normal distributions with mean $\mu _{k,l,{stoch}}$ and variance $\sigma ^2_{k,l}$. In the deterministic case, the observation consists of these fields with Fourier coefficients of prescribed magnitudes $\mu _{k,l, {det}}$. The prediction is subsequently nudged towards the observation by correcting the predicted magnitudes of the Fourier coefficients. The weight of the model, often referred to as the ‘gain’, is determined by $\alpha _j\Delta t / \tau _{k, l}$ for each $k, l$ separately.

3.2. Heat transport correction

The heat transport in the turbulent flow is described by the Nusselt number and is considered the key response of the system to the imposed Rayleigh number (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009). The definition of the Nusselt number that we adopt here is

(3.8)\begin{equation} Nu = 1 + \sqrt{Pr Ra}\langle v T\rangle_\varOmega, \end{equation}

which is well suited for use on coarse computational grids. An alternative definition of $Nu$ involves a gradient of temperature, which is more sensitive to coarse-graining. In (3.8), $\varOmega$ denotes the domain with area $|\varOmega |$ and $\langle \boldsymbol {\cdot }\rangle _\varOmega$ denotes the domain average. It is clear from definition (3.8) that $vT$ needs to be modelled accurately to recover skilful predictions of the heat flux. To achieve this, we propose a constraint to be used in conjunction with the model described in § 3.1.

The volume average in (3.8) comprises averages of the heat flux along horizontal cross-sections of the domain. For a fixed vertical coordinate $y_l$, we denote the heat flux along this cross-section by $\langle v T \rangle _l$. Along this cross-section, we indicate the Fourier coefficients of the velocity and temperature with a hat symbol $\hat {\cdot}$ and observe that

(3.9)\begin{equation} (\widehat{vT})_0 = \sum_{k}\hat{v}^*_k\hat{T}_k, \end{equation}

where the subscript $k$ signifies the $k^\mathrm {th}$ Fourier coefficient. The subscript $0$ indicates that we consider the zeroth mode of the Fourier series, which by definition equals the value of $\langle v T\rangle _l$. We assume that a mean heat flux along the horizontal cross-section is known from the reference high-resolution data and denote this value by $F_l$. Subsequently, the heat flux along the horizontal cross-section in a coarse numerical simulation is corrected by minimizing the error $\|F_l - \sum _k \hat {v}_k^* \hat {T}_k \|^2$ with respect to $T$. Here, we minimize the error by varying the phases of the Fourier coefficients $\hat {T}_k$. We alter the temperature instead of the vertical velocity so that the velocity field remains divergence-free. Adapting the phases only ensures that the spectrum of the temperature along the horizontal cross-section is invariant under the heat transport correction. In total, the heat flux correction is an extension of the nudging procedure (3.6). It enables a correction of the temperature field solely based on a statistic of the reference solution, rather than on a dynamic equation. In doing so, the dependence between the vertical velocity and the temperature is taken into account in the nudging procedure. Thus, applying the heat flux correction ensures an improved average Nusselt number estimate and is therefore expected to improve the accuracy of the numerical solutions.

The error $\|F_l - \sum _k \hat {v}_k^* \hat {T}_k \|^2$ is minimized using a gradient descent algorithm. We note that the correction may in principle yield an arbitrarily good approximation of the reference heat flux, but this is not guaranteed to produce physically relevant results. Instead of aiming for an exact agreement of the mean heat flux, we apply the gradient descent algorithm until the heat flux in the horizontal cross-section is within a 10 % margin of the reference value. This serves to demonstrate the added value of the correction. Preliminary tests have shown that this already improves the heat flux significantly without qualitatively altering the temperature field. In the next section, coarse numerical simulations are performed both with and without the heat flux correction. The optimization of this procedure is beyond the scope of this paper.

4.1. Qualitative model performance

A qualitative comparison of the model configurations M0-7 is given in figures 3–6 by means of instantaneous snapshots of the obtained statistically stationary states. In these figures, the snapshots of the DNS and of M0 are the same as depicted in figure 1. A comparison of the temperature fields in statistically steady states is provided in figure 3. Here, we observe that the configurations M1–4 do not lead to significant qualitative changes in the temperature field when compared with the no-model configuration M0. In these configurations, the temperature is not explicitly forced and suffers from artificial dissipation inherent to the coarsening. The model configurations M5–7, in which the temperature is forced directly, display more pronounced small-scale features. At the same time, the large-scale circulation pattern is still visible in these results. In addition, from the results of M6 and M7, we conclude that applying the heat flux correction does not lead to qualitatively different temperature fields.

Figure 3. Temperature fields in statistically steady states. Shown are the reference solution and the results obtained with coarse numerical simulations M0–7. The colour scheme is the same as used for the temperature fields shown in figure 1. (a) DNS, (b) M0, (c) M1, (d) M2, (e) M3, ( f) M4, (g) M5, (h) M6, (i) M7.

The pressure fields of the corresponding solutions are shown in figure 4. Here, we recall that any detail observed in the DNS pressure field is lost in the coarse no-model result M0. No improvements are observed in the pressure field when only the temperature is explicitly forced, as is done in M5. The remaining model configurations all yield a distinct qualitative improvement in the pressure fields. In particular, only applying a large-scale velocity correction already qualitatively changes the pressure field. This is observed for the deterministic and the stochastic forcing, given by M1 and M3, respectively. The addition of forcing the velocity at small scales or simultaneously forcing the temperature does not yield additional significant changes. A noticeable difference exists between the deterministic and stochastic methods. As becomes clear from M3–4, the random forcing leads to a fragmentation of the coherent structures in the pressure field.

Figure 4. Pressure fields in statistically steady states. Shown are the reference solution and the results obtained with coarse numerical simulations M0–7. The colour scheme is the same as used for the pressure fields shown in figure 1. (a) DNS, (b) M0, (c) M1, (d) M2, (e) M3, ( f) M4, (g) M5, (h) M6, (i) M7.

The horizontal velocity fields and vertical velocity fields are provided in figures 5 and 6, respectively. We observe that all coarse numerical solutions display agreement with the DNS in terms of large-scale coherent structures. Nonetheless, artificial dissipation leads to an underestimate of the velocity magnitude in cases M0 and M5. This suggests that only forcing the temperature is not sufficient for accurately reproducing the velocity fields. The other cases indicate that explicitly forcing the velocity leads to accurate velocity magnitudes.

Figure 5. Horizontal velocity fields in statistically steady states. Shown are the reference solution and the results obtained with coarse numerical simulations M0–7. The colour scheme is the same as used for the horizontal velocity fields shown in figure 1. (a) DNS, (b) M0, (c) M1, (d) M2, (e) M3, ( f) M4, (g) M5, (h) M6, (i) M7.

Figure 6. Vertical velocity fields in statistically steady states. Shown are the reference solution and the results obtained with coarse numerical simulations M0–7. The colour scheme is the same as used for the vertical velocity fields fields shown in figure 1. (a) DNS, (b) M0, (c) M1, (d) M2, (e) M3, ( f) M4, (g) M5, (h) M6, (i) M7.

4.2. Energy spectra

We now establish that the model proposed in § 3.1 improves the average energy spectra of the forced variables. The average energy spectra of the velocity components and the temperature are shown in figure 7, displaying the spectra along a horizontal cross-section near the wall and in the centre of the domain. Both near the wall and in the centre of the domain, respectively shown in figure 7(a–c) and figure 7(d–f), the no-model M0 results exhibit significant differences compared with the filtered DNS. The measured energy levels of the velocity are too low with M0 at all resolved scales. In contrast, a significant discrepancy in the temperature spectra is observed only for wavenumbers larger than 10.

Figure 7. Time-averaged energy spectra measured along horizontal cross-sections of the domain for (a,d) the horizontal velocity, (b,e) vertical velocity and (c, f) temperature. The cross-sections are taken near (a–c) the bottom wall and (d–f) the core of the domain. The cross-sections are taken at $y=8.5\times 10^{-4}, y=5.5\times 10^{-1}$ for the horizontal velocity, and at $y=5.0\times 10^{-4}, y=5.0\times 10^{-1}$ for the vertical velocity and the temperature. The solid lines show the average spectra of the filtered DNS, the model results are displayed using the symbols in table 1.

The discrepancies in the spectra of M0 and the reference are attributed to artificial dissipation caused by the coarsening. In particular, the numerical dissipation affects both the velocity and the temperature spectra at higher wavenumbers. Through the nonlinear interactions in the momentum equation, the velocity is adversely affected at all wavenumbers. This is further corroborated by the results of M2 and M4, where all available length scales are forced only for the velocity. In the core of the domain, where the coarsening is strongest, these results display accurate velocity spectra but yield no improvement in the temperature spectra, suggesting that the temperature still suffers from artificial viscosity in these cases. Apparently, the improvements in the velocity spectra influence the prediction of the temperature only to a small degree.

The large-scale velocity forcing applied in M1 and M3 yields improved velocity energy levels at low wavenumbers. However, the improvement gradually vanishes at higher wavenumbers. These configurations exhibit no improvement in the temperature spectra. The cases M2 and M4 lead to an improved agreement on the velocity spectra at all wavenumbers in the centre of the domain, establishing the spectrum-reconstructing property of the model described in § 3.1. Nonetheless, all models underestimate the large-scale energy in the centre of the domain. At these scales, the measured correlation time $\tau$ is large and therefore the model contribution is limited.

Near the wall, the horizontal velocity is accurately represented at all wavenumbers despite the fact that the energy of the vertical velocity deviates from the reference for wavenumbers larger than 15. The temperature spectra for M2 and M4 near the wall show good agreement with the reference. Comparing this with the results of M1 and M3 indicates that the prediction of near-wall temperature is improved by the forcing of small-scale velocity despite no explicit forcing being applied to the temperature. No improvement is observed for these cases in the centre of the domain, which we attribute to artificial dissipation.

The velocity spectra show no significant change when only the temperature is explicitly forced, as observed from the results of M5. This case produces an accurate temperature spectrum in the core of the domain and yields an improved spectrum near the wall. Additionally forcing the velocity significantly improves the velocity spectra, as is observed for cases M6 and M7. Here, we observe good agreement for the velocity and the temperature across all length scales in the centre of the domain. In particular, a definite improvement is observed when comparing the temperature spectrum to those of M1–4. Near the wall, the horizontal velocity and the temperature are both captured accurately, while the vertical velocity still deviates for wavenumbers larger than 15. The similarity between the spectra obtained for M6 and M7 indicates that the heat flux correction described in § 3.2 does indeed not lead to significant changes in the spectra.

4.3. Flow statistics

The mean temperature and mean heat flux are displayed in figure 8 as a function of the wall-normal distance. All models except M5 efficiently mitigate the small mean temperature discrepancy between M0 and the reference.

Figure 8. Comparison of the (a) time-averaged temperature and (b) time-averaged heat flux measured along horizontal cross-sections of the domain and displayed as a function of the wall-normal distance. The solid line shows the mean values of the filtered DNS, the model results are displayed using the symbols in table 1.

The mean heat flux of the no-model M0 case is consistently too low, which is a direct result of underestimating the vertical velocity. Applying the large-scale velocity forcing as done in cases M1 and M3 yields an improved heat flux. In particular, the measured heat flux near the wall shows good agreement with the reference. The mean heat flux is consistently overestimated when the velocity is forced at all wavenumbers, which is the case for M2, M4 and M6. Comparison of the results of M6 with M7 establishes that the heat flux correction described in § 3.2 ensures a better prediction of the mean heat flux. Finally, only imposing the temperature spectrum deteriorates the measured heat flux, as shown by the results of M5.

These observations in combination with the energy spectra of the previous subsection expose the simplifying model assumptions discussed in § 3.1. Despite accurate energy spectra of all variables, the M6 model does not yield an accurate heat flux. This indeed suggests that the energy spectra alone do not provide sufficiently strict modelling criteria for obtaining accurate coarse-grid numerical simulations, and instead benefit from additional cross-variable constraints such as the imposed heat flux.

The r.m.s. of the velocity components are shown in figure 9(a,b) as a function of the wall-normal distance. A strong reduction of the turbulent intensity of the velocity is observed for the no-model M0 results. Similar to previous observations for case M5, only forcing the temperature does not lead to improvements in the r.m.s. of the velocity. All other model configurations lead to a comparable improvement in the r.m.s. of the horizontal velocity. A slight difference between the stochastically forced and deterministically forced solutions may be distinguished in the r.m.s. profiles of the horizontal velocity, visible in the results of M3–4. Comparable results are observed for the r.m.s. of the vertical velocity, where all models except M0 and M5 display good agreement with the reference.

Figure 9. Root mean square (r.m.s.) of the (a) horizontal velocity, (b) vertical velocity and (c) temperature, measured along horizontal cross-sections of the domain and displayed as a function of the wall-normal distance. The solid line shows the r.m.s. values of the filtered DNS, the model results are displayed using the symbols in table 1.

Figure 10. Comparison of the rolling mean of the kinetic energy (KE) over time. The solid line shows the KE of the filtered DNS, and the model results are displayed using the symbols in table 1.

The average temperature fluctuations are shown in figure 9(c). We observe that all model configurations except M0 and M5 predict the wall-normal distance of the peak of the fluctuations accurately. However, the model overestimates the maximal predicted r.m.s. by $7.5\,\%$ to $18\,\%$.

4.4. Total kinetic energy and heat flux

A comparison of the rolling mean of the total kinetic energy (KE) is shown in figure 10. The improvement obtained by M1–4, M6 and M7 is evident. The coarse simulations are initialized from a snapshot of the filtered DNS in the reference statistically stationary state. Due to coarsening effects, this statistically stationary state cannot be sustained in coarse simulations without including a model term. Numerical dissipation and other discretization effects lead to deviations from the filtered DNS on coarse grids, despite adopting the same physical parameters as in the high-resolution simulation. The forcing model aims to steer the coarse numerical simulation towards the reference statistically stationary state by applying a correction at each time step. Nonetheless, this procedure is approximate and does not necessarily yield the same total kinetic energy and Nusselt number as observed in the reference. Therefore, transient behaviour is observed initially until the coarse simulations settle at a different statistically steady state. At $t=400$, the mean of the KE for M0 is approximately $31\,\%$ of the reference KE. Only forcing the temperature, shown by M5, deteriorates the total energy and yields roughly $27\,\%$ of the reference value. The other models contain between $72\,\%$ and $77\,\%$ of the reference value. It is reasonable to assume that this discrepancy is predominantly caused by the model underestimating the energy in large scales in the centre of the domain, as was discussed in § 4.2.

A quantification of the total heat flux in the domain is provided by comparing the time-averaged Nusselt number, shown in figure 11. Note that the reference value ${Nu=95}$ is shown with $5\,\%$ error margins. The no-model coarse-grid simulation leads to an underestimated heat flux resulting from the reduced velocity magnitude induced by artificial dissipation. The temperature forcing in case M5 was previously shown to not yield any improvements in the mean temperature or the velocity and does therefore not improve the Nusselt number estimate. A correction of the large-scale velocity features in configurations M1 and M3 leads to a very accurate Nusselt number estimate. Nonetheless, an accurate description of the velocity does not guarantee an accurate heat flux. This is underpinned by the results of M2, M4 and M6, which all exhibit an accurate representation of large and small velocity features, but consistently overestimate the Nusselt number. Finally, we observe that this adverse effect is efficiently mitigated by the heat flux correction, as becomes evident from the resulting Nusselt number estimate of M7.

Figure 11. Comparison of the rolling mean of the Nusselt number over time. The solid line at $Nu=95$ shows the theoretically predicted value, with $5\,\%$ error margins given by the dashed lines. The model results are displayed using the symbols in table 1.

4.5. Robustness under variations in forcing strength

The robustness of flow predictions with respect to the forcing strength is addressed next. A comparison of various forcing strengths is carried out by multiplying all measured $\tau _{k,l}$ by $0.5$, increasing the forcing strength, or by 2, 4 and 8, decreasing the forcing strength. These scalings are applied to the model configurations M0–7 described in table 1. Figures 12 and 13 show the energy spectra in the bulk of the domain and the r.m.s. of the function of the wall-normal distance, respectively, obtained when multiplying the measured $\tau _{k,l}$ by 8. The results are only shown for this scaling since these differ the most from the results obtained with the unscaled forcing strength.

Figure 12. Time-averaged energy spectra measured along horizontal cross-sections of the domain for the (a) horizontal velocity, (b) vertical velocity and (c) temperature. The cross-sections are taken near the core of the domain, at $y=5.5\times 10^{-1}$ for the horizontal velocity, and at $y=5.0\times 10^{-1}$ for the vertical velocity and the temperature. The model results are obtained by multiplying the measured correlation times by 8. The solid lines show the average spectra of the filtered DNS, the model results are displayed using the symbols in table 1.

Figure 13. Root mean square (r.m.s.) of the (a) horizontal velocity, (b) vertical velocity and (c) temperature, measured along horizontal cross-sections of the domain and displayed as a function of the wall-normal distance. The model results are obtained by multiplying the measured correlation times by 8. The solid line shows the r.m.s. values of the filtered DNS, the model results are displayed using the symbols in table 1.

Increasing the forcing strength leads to numerical solutions that closely follow the reference kinetic energy spectra. As a result, the r.m.s. values as a function of the wall-normal distance also improve considerably with respect to the no-model coarse numerical simulation. However, no significant differences have been observed when compared with the results obtained with the original (unscaled) values of $\tau _{k,l}$. Decreasing the forcing strength leads to coarse numerical simulations that follow the reference energy spectra less closely, but nonetheless improve upon the no-model coarse simulation. This is displayed in figure 12. This result is also reflected in the r.m.s. as illustrated in figure 13.

Changing the forcing strength has a noticeable effect on the total kinetic energy and the measured Nusselt number. The time-averages of the total kinetic energy and Nusselt number are depicted in figure 14 as a function of the scaling of the forcing strength. All model configurations except M5, where only the temperature is explicitly forced, show a decrease of the kinetic energy and of the Nusselt number as the forcing strength is decreased.

Figure 14. Comparison of the final rolling mean values of the (a) total kinetic energy (KE) and (b) Nusselt number for for different scalings of the correlation times $\tau _{k,l}$ used in the model. The model results are displayed using the symbols in table 1.

4.6. Generalization to other Rayleigh numbers

We now turn our attention to the model performance at different Rayleigh numbers. Increasing the Rayleigh number in three-dimensional Rayleigh–Bénard convection affects the velocity spectra and the temperature spectra along horizontal cross-sections. De, Eswaran & Mishra (Reference De, Eswaran and Mishra2017) considered a range of Rayleigh numbers between $7\times 10^4$ and $2\times 10^6$ at $Pr=0.71$ in a large aspect ratio box. An increase of high-wavenumber excitation is reported for kinetic energy, temperature fluctuations, heat flux and kinetic energy dissipation, and is observed near the top and bottom plates and in the bulk of the domain. The maximum energy dissipation is found to shift towards higher wavenumbers as the Rayleigh number is increased. This might suggest that a coarse simulation at larger Rayleigh numbers benefits from increased forcing at larger wavenumbers.

Figure 15. Measured Nusselt numbers as a function of the Rayleigh number. The solid red line shows reference Nusselt number predictions from the literature, with 5 % error margins given by the dashed lines. The no-model and large-scale forcing results are obtained with configurations M0 and M1, respectively. The reference scaling exponents are $0.285$ for $10^9\leq {Ra}\leq 10^{13}$ and $0.35$ for ${Ra}\geq 10^{13}$.

When increasing the Rayleigh number, the mean velocity and temperature profiles in two-dimensional Rayleigh–Bénard convection behave similarly when normalized by the friction velocity and corresponding characteristic temperature, respectively (Zhu et al. Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018). In the range between $Ra=10^{11}$ and $Ra=10^{14}$, the mean profiles displayed the same behaviour in the viscous sublayer near the wall at each simulated Rayleigh number. This is followed by a log layer for the velocity and a flat region for the temperature. These results indicate that the mean profiles do undergo quantitative changes but no qualitative changes in this range of Rayleigh numbers. Therefore, a model calibrated at a particular Rayleigh number might still have merit when different Rayleigh numbers are considered in coarse numerical simulations.

Here, we assess the Nusselt number estimates at different Rayleigh numbers using the forcing calibrated at $Ra=10^{10}$. The Nusselt number estimates in coarse simulations without a model (M0) and with deterministic large-scale momentum forcing (M1) are compared with reference values reported by Johnston & Doering (Reference Johnston and Doering2009) and Zhu et al. (Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018). The reference scaling exponents for the Nusselt number as a function of the Rayleigh number are $0.285$ for $10^9\leq {Ra}\leq 10^{13}$ and $0.35$ for ${Ra}\geq 10^{13}$. We consider only the model configuration M1 in this comparison, since this configuration was already found to produce an accurate Nusselt number at $Ra=10^{10}$ without using explicit knowledge of the heat flux. The results are summarized in figure 15, showing the measured Nusselt number as a function of the Rayleigh number. These results show that the large-scale momentum forcing yields accurate Nusselt number estimates across several decades of Rayleigh numbers without using high-resolution data for these parameter regimes. At the transition to the so-called ultimate regime, observed at $Ra=10^{13}$ (Zhu et al. Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018), the Nusselt number estimates of the coarse numerical simulations deteriorate somewhat. The estimates lose accuracy for much larger Rayleigh numbers compared with the reference $Ra=10^{10}$. This is likely a result of the gradually accumulating quantitative flow differences as the Rayleigh number increases. Extending the range of accurate predictions provides a challenge for future work.