Data-driven spectral turbulence modelling for Rayleigh–Bénard convection

Abstract

A data-driven turbulence model for coarse-grained numerical simulations of two-dimensional Rayleigh–Bénard convection is proposed. The model starts from high-fidelity data and is based on adjusting the Fourier coefficients of the numerical solution, with the aim of accurately reproducing the kinetic energy spectra as seen in the high-fidelity reference findings. No assumptions about the underlying partial differential equation or numerical discretization are used in the formulation of the model. We also develop a constraint on the heat flux to guarantee accurate Nusselt number estimates on coarse computational grids and high Rayleigh numbers. Model performance is assessed in coarse numerical simulations at $Ra=10^{10}$. We focus on key features including kinetic energy spectra, wall-normal flow statistics and global flow statistics. The method of data-driven modelling of flow dynamics is found to reproduce the reference kinetic energy spectra well across all scales and yields good results for flow statistics and average heat transfer, leading to computationally cheap surrogate models. Large-scale forcing extracted from the high-fidelity simulation leads to accurate Nusselt number predictions across two decades of Rayleigh numbers, centred around the targeted reference at $Ra=10^{10}$.

1. Introduction

Turbulent flows are characterized by the distribution of kinetic energy over a vast range of scales. The nonlinearity in the Navier–Stokes equations ensures that large and small eddies interact with each other, resulting in a wide range of dynamic flow features (Pope & Pope 2000). This process transfers kinetic energy from the large energy-containing scales towards smaller scales, until the kinetic energy is finally dissipated by viscosity at the smallest scales. The energy cascade towards the smallest scales yields a significant challenge in computational fluid dynamics in the turbulent regime (Geurts 2003; Sagaut 2006). To accurately simulate the flow, the fluid dynamical model should resolve the scales of turbulence down to the Kolmogorov length scale. A direct approach would then require very fine computational grids which is often intractable even with modern-day computing resources. A common way to alleviate the large computational requirements is by reducing the numerical resolution at which an approximate solution to the flow is obtained. To compensate for the lack of refinement of the computational approach, a model term is subsequently added to the governing equations to represent the influence of unresolved dynamics on the resolved components of the flow (Geurts & Holm 2002; Piomelli, Rouhi & Geurts 2015; Rouhi, Piomelli & Geurts 2016; Geurts 2022).

In this paper, we describe how prior knowledge of flow statistics obtained from an offline fully resolved simulation may be incorporated to construct an online high-fidelity model for coarse numerical simulations. The proposed reduced-order model acts on the numerical solution in spectral space, employing techniques from time series modelling and data assimilation. This model is designed to yield accurate kinetic energy spectra, despite the rather coarse flow representation. We demonstrate the capabilities of this data-driven approach for two-dimensional Rayleigh–Bénard convection at high Rayleigh number.

Rayleigh–Bénard convection is a fundamental problem in fluid dynamics, describing buoyancy-driven flows heated from below and cooled from above (Kadanoff 2001; Ahlers, Grossmann & Lohse 2009; Kunnen, Geurts & Clercx 2009; Kooij et al. 2018). In particular, thermal convection is meaningful in geophysical processes, such as in describing convective processes in the atmosphere (Hartmann, Moy & Fu 2001) or the ocean (Marshall & Schott 1999). The large range of scales present in turbulence is also further affected by buoyancy effects. For example, a common phenomenon observed in Rayleigh–Bénard convection is the formation of spatially coherent structures in large-scale circulation (Ahlers et al. 2009) and, in larger spatial domains, the formation of thermal superstructures (Stevens et al. 2018). However, a thin boundary layer exists near either of the walls which becomes turbulent and increasingly thinner with growing temperature differences between the walls, i.e. growing Rayleigh number (Kraichnan 1962; Zhu et al. 2018). Properly resolving the boundary layers requires large computational grids and poses a challenge even in two-dimensional Rayleigh–Bénard convection (Zhu et al. 2018). This stresses the conundrum of computational fluid dynamics, where one strives to find a balance between simulating flows at modest computational costs without adversely affecting the prediction of flow statistics.

Simulating flows at modest computational costs while retaining a high level of accuracy is the aim of large-eddy simulation (LES) (Geurts 2003; Sagaut 2006). Instead of fully resolving all length scales of the flow, a computationally less intensive approximation is found by coarsening the flow description and simultaneously including a subgrid model to accommodate the loss of explicit finer details in the dynamics. The coarsening is accomplished by spatial filtering, which, through the specification of the filter width, establishes the required level of refinement that should be included in the numerical simulations. The subgrid model then approximates the effect the unresolved dynamics has on the resolved scales, and serves as a closure for the filtered equations. This approximation depends on both the adopted filter (Geurts & Holm 2003) and selected closure model, as well as the choice of discretization and level of coarsening (Langford & Moser 1999; van der Bos, van der Vegt & Geurts 2007; Beck, Flad & Munz 2019).

With the increase of computational resources, direct numerical simulation (DNS) of turbulent flows is achievable to an ever-increasing extent and may serve to generate data from which LES models could be derived. This data-driven LES has been an active topic of research in recent years. For example, the decomposition of unresolved dynamics into fixed global basis functions and corresponding time series yields an efficient approach for which only the latter needs to be modelled. Frederiksen & Kepert (2006) employed DNS data of the barotropic vorticity equation to model the time series of spherical harmonics as stochastic processes with memory effects, leading to accurate kinetic energy spectra in coarse-grid simulations. Using proper orthogonal decomposition (POD), Ephrati et al. (2022b) showed that applying corrections to coarse-grid numerical simulations may lead to significant error reduction. Machine-learning methods have also been successfully employed to find subgrid models (Beck et al. 2019), reporting improved results compared to traditional eddy-viscosity models. Examples include using artificial neural networks in two-dimensional decaying turbulence (Maulik et al. 2019) and convolutional neural networks in three-dimensional homogeneous isotropic turbulence (Kurz, Offenhäuser & Beck 2023), yielding improved energy spectra and turbulent fluctuation distributions. Similar methods have also been applied to two-dimensional Rayleigh–Bénard convection. A reduced-order model is developed by Pandey et al. (2022) employing machine learning, where a combined convolutional autoencoder–recurrent neural network (Vlachas et al. 2020) is used to predict the local turbulent heat flux. Accurate probability density functions are obtained for the local heat flux, while a dimensionality reduction to 0.2 % of the original size is achieved for the turbulence data. In a similar vein, Heyder & Schumacher (2021) study moist Rayleigh–Bénard convection. The dimensionality of the system is reduced by applying POD to the high-fidelity data, after which the POD coefficients are predicted using echo state networks. Here a good agreement of low-order flow statistics is observed in the reduced-order model.

Incorporating data into numerical models to improve flow predictions is well established in geophysical fluid dynamics, where data assimilation has been successfully employed for several decades. The aim is to improve forecasting by minimizing the differences between observed and modelled values while accounting for uncertainties (Ghil & Malanotte-Rizzoli 1991; Daley 1992). In particular, continuous data assimilation (CDA) aims to nudge the model solution towards an observed reference by means of a feedback control term acting as external forcing (Azouani & Titi 2013; Azouani, Olson & Titi 2014). This concept is also extended to linear stochastic differential equations, arising as the continuous-time limit of the three-dimensional variational (3DVAR) data assimilation algorithm (Blömker et al. 2013), which has been shown to successfully steer coarse-grained numerical solutions of the two-dimensional Navier–Stokes equations towards an observed reference solution. Additionally, the convergence of coarse numerical solutions augmented by CDA to an observed reference solution has been proven (Farhat, Jolly & Titi 2015) and shown numerically for two-dimensional Rayleigh–Bénard convection (Altaf et al. 2017).

The purpose of this paper is to combine ideas from data assimilation with large-eddy simulation. In particular, we derive a model term based on statistical quantities from a reference high-resolution simulation and use this as a stand-alone model for coarse numerical simulations. Our proposed method incorporates Ornstein–Uhlenbeck processes in the evolution of the Fourier coefficients of the numerical solution, steering the solution towards an a priori measured statistically steady state. Only three parameters need to be defined for each Fourier mode, outlining the simplicity of the model. The parameters are inferred from data, do not depend on the adopted spatial or temporal discretization and are defined such that the reference energy spectrum is closely reproduced in the coarse-grid simulations. The resulting prediction-correction scheme is of the form of the diagonal Fourier domain Kalman filter (Harlim & Majda 2008; Majda & Harlim 2012) with a fixed prescribed gain. This identification enables future research that combines LES and data assimilation. The same approach has been applied in a recent study of coarse-grid modelling of the two-dimensional Euler equations on the sphere (Ephrati et al. 2023a), where a decomposition of the vorticity field into spherical harmonic basis modes was employed in the coarse-grid model.

The two main results reported in this paper are the following. First, the high-resolution data at $Ra=10^{10}$ can be used successfully to define a forcing acting on large scales of coarse numerical simulations. By improving the energy content in these forced scales in coarsened simulations, an accurate estimate of the Nusselt number is obtained. The forcing calibrated at $Ra=10^{10}$ is found to yield accurate Nusselt number estimates across two decades of Rayleigh numbers, from $Ra=10^{9}$ to $Ra=10^{11}$, without the need to re-compute a high-fidelity simulation for each Rayleigh number of interest. Second, combining forcing across all resolvable scales with prior knowledge of the heat flux in the domain leads to an ‘offline/online’ approach through which an accurate and effective coarsened model for online real-time predictions can be formulated. The model that yields such accurate coarsening is derived from the high-fidelity offline simulation. The result is a computationally cheap surrogate model for Rayleigh–Bénard convection, accurately representing the complex behaviour up to the smallest resolvable scales on the coarse grid. This is demonstrated here at $Ra=10^{10}$.

The paper is structured as follows. The governing equations and adopted discretization are described in § 2. The data-driven model is introduced in § 3, detailing the forcing in § 3.1 and complementary heat flux correction in § 3.2. The model performance for a variety of model configurations is assessed in § 4. Conclusions are presented in § 5.

2. Governing equations and numerical methods

In this section, we introduce the governing equations and the simulated case in § 2.1, the employed numerical discretization in § 2.2 and illustrate the effects of severe coarsening, without any modelling correction, on the quality of the solution in § 2.3.

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. 2018). 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 2013).

2.2. Numerical methods

The adopted spatial discretization is an energy-conserving finite difference method (Vreman 2014) and is parallelized as by Cifani, Kuerten & Geurts (2018). 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. 2015). 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. (2015). 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 1991; van der Poel et al. 2015; Cifani et al. 2018). 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 1985). 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 1979). 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. 2018). 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 2005). 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. Spectrum-preserving forcing

In this section, we describe a data-driven forcing to augment coarsened numerical simulations of statistically steady states. The high-resolution and low-resolution snapshots presented in the previous section hinted at the need of introducing explicit forcing to more accurately represent average flow patterns on coarse computational grids. Here, we propose a model to reproduce the kinetic energy spectra in coarse numerical simulations.

The model parameters are extracted from the reference data, obtained from a sequence of 8040 solution snapshots each separated by $0.05$ time units of a DNS performed at a $4096\times 2048$ computational grid. With these numerical data, we achieve sufficiently many snapshots to reliably recover statistical properties of the flow, which is a pre-requisite for our model development. The Fourier components of horizontal cross-sections of the solution are computed for each snapshot, yielding time series for each streamwise wavenumber at each $y$-coordinate. The magnitudes of these complex time series yield mean values, variances and correlation times that are used as model parameters. We next present the main steps in the construction of this model.

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. 2023b), although the common choice is to let it be a Gaussian process with a variance depending on the time step size (Higham 2001). 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. 2013). 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 2005) 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 2008; Majda & Harlim 2012). 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. 2009). 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. Model performance

In this section, we apply the model in eight different configurations to numerical simulations on the coarse grid. The configurations are listed in table 1 and differ in the variable that is being forced, the wavelengths at which the forcing is applied and whether the forcing is deterministic or stochastic. These configurations will be referred to as M0-7, inspired by the nomenclature used in the comparison of LES models by Vreman, Geurts & Kuerten (1997). Here, the wavenumbers at which the forcing is applied are chosen as $l\leq 5$ and $l\leq 32$. The former implies that the model only explicitly acts on the large scales of motion and the latter implies that all resolved scales are directly affected by the model. This set of configurations is chosen to distinguish between the effects of large-scale forcing and small-scale forcing, deterministic forcing and stochastic forcing, and the choice of the forced variable. The model simulations are run with a time step size $\Delta t = 0.02$. Recall that the filtered DNS provides the reference solution.

Table 1. Model configurations used in the coarse numerical simulations.

Figure 2 shows the correlation times of the Fourier coefficients $\tau _{k,l}$ for each pair $k, l$ for the horizontal velocity, the vertical velocity and the temperature, measured from the high-fidelity data. As expected, large streamwise wavenumbers, i.e. structures with small length scales, correspond to small correlation times. Conversely, small wavenumbers, i.e. large length scale structures, correspond to large correlation times. This confirms that the model contribution to the dynamics of the Fourier coefficients at small scales is stronger than at large scales. The horizontal velocity and the temperature generally yield smaller correlation times in the bulk of the domain compared to the region near the walls, whereas the opposite is observed for the vertical velocity. There is a notable difference at the zero-wavenumber correlation time for the vertical velocity, which is due to the fact that the average vertical velocity is always zero. Finally, we note that $\tau _{k,l}$ is observed to be nowhere below $\alpha _j\Delta t$, thus implying that the magnitudes of the Fourier coefficients do not depend solely on the model in any case, but also take contributions from the discretized equations directly.

Figure 2. Correlation times $\tau _{k,l}$ of the Fourier coefficients for the (a) horizontal velocity, (b) vertical velocity and (c) temperature measured from the high-fidelity data.

We first provide in § 4.1 an impression of the qualitative improvements obtained when applying the model. In the ensuing subsections, a detailed quantitative comparison is carried out. Several quantities will be compared with the filtered DNS data to gain insight into the quality of the model. In § 4.2, we first verify that the model approximates the average energy spectra of the filtered DNS by comparing the spectra of the velocity and the temperature near the wall and in the core of the domain. In § 4.3, the mean temperature, the mean heat flux and the root-mean-square deviation (r.m.s.) are measured as a function of wall-normal distance and compared with the reference. Finally, global flow statistics such as the total kinetic energy and the Nusselt number are examined in § 4.4.

The r.m.s., mean temperature and mean heat flux rely on averages along horizontal cross-sections of the domain. For a fixed value of $y$, we adopt the following definition:

(4.1)\begin{equation} \mathrm{r.m.s.}\ (f, y, t) = \left[\frac{1}{|A|}\int_A \left(f(x, y, t) - \langle f(x, y, t)\rangle_A \right)^2\,\text{d}A \right]^{1/2}, \end{equation}

where $\langle \boldsymbol {\cdot }\rangle _A$ denotes the average over the horizontal cross-section with length $|A|$ and $f$ is the field of interest. The mean temperature and heat flux are computed as the mean $\langle \boldsymbol {\cdot }\rangle _A$ of the corresponding fields. The global kinetic energy will be computed as

(4.2)\begin{equation} \mathrm{KE} = \int_\varOmega \tfrac{1}{2}(u^2 + v^2)\,\text{d}\varOmega \end{equation}

and the Nusselt number follows from (3.8).

Our interest lies in the time average of the aforementioned quantities. The quality of coarse-grid models is therefore measured by comparing averaged quantities rather than instantaneous quantities (Vreman et al. 1997; Langford & Moser 1999). The energy spectra, r.m.s. values and mean temperature, and heat flux will be measured after the coarse-grid numerical simulations have reached a statistically steady state. The global quantities of interest are illustrated using a rolling average over time.

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 (2017) 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. 2018). 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 (2009) and Zhu et al. (2018). 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. 2018), 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.

5. Conclusions and outlook

In this paper, we have proposed a data-driven model for coarse numerical fluid simulations and assessed its performance when applied to two-dimensional Rayleigh–Bénard convection. Statistical information of Fourier coefficients of a reference direct numerical simulation was used to infer model parameters, which constituted a forcing term for reproducing the reference energy spectra. The model parameters are defined such that small scales of motion are corrected strongly. Various model configurations were applied to gain insight into the model performance, generally leading to improved results compared with using no model.

Applying the model at all wavelengths resulted in significant improvement of the spectra both near the walls and near the centre of the domain, which established that the model had its desired effect on the numerical solution. Additionally, the application of the model was found to yield improved estimates of flow statistics. In particular, the average turbulent fluctuations and average temperature improved significantly compared with the no-model case. The total kinetic energy was found to improve upon using the model but highlighted that modelling the effects of unresolved dynamics on the large-scale flow features as independent processes might be too restrictive. The measured total heat flux was accurately captured for several model configurations, although accurately reconstructing energy spectra was shown not to be a sufficient criterion for this purpose. The latter problem was efficiently alleviated by including a constraint on the average heat flux in the model, leading to a computationally cheap data-driven surrogate model for the highly complex flow dynamics. Applying large-scale momentum forcing in the coarse numerical simulations yielded an accurate Nusselt number estimate without requiring explicit knowledge of the heat flux. Furthermore, accurate Nusselt number estimates were obtained across two decades of Rayleigh numbers using only the forcing calibrated at $Ra=10^{10}$.

Future work will be dedicated to expanding the proposed model by consulting Kalman filtering theory. Specifically, the interactions between the Fourier modes can be explicitly represented by including covariance estimates in the model. This would additionally serve to verify at which frequencies the Fourier modes evolve independently, which is expected to result in a better understanding of the modelling of small-scale flow features. Alternatively, spatially coherent structures can be included in the model by applying proper orthogonal decomposition to the reference data, as demonstrated by Ephrati, Cifani & Geurts (2022a). Although no assumptions are made about the numerical method or adopted coarse resolution in the formulation of the model, further numerical experiments adopting a different resolution or discretization may be carried out to assess the robustness of the model. Finally, we note that the presented method might be extended to three-dimensional flows. This could be achieved by extending the method presented in this paper to also cover the reconstruction of two-dimensional energy spectra. For simple domains, covered with a structured tensor grid, this could be achieved by treating each horizontal (or constant-index) cross-section of the three-dimensional domain separately. Alternatively, one can employ proper orthogonal decomposition and recover the POD spectrum as done by Ephrati et al. (2022a) and readily extend this from two dimensions to three dimensions.