# Continuous and discrete data assimilation with noisy observations for the Rayleigh-Bénard convection: a computational study

Mohamad Abed El Rahman Hammoud, Olivier Le Maître, Edriss S Titi, Ibrahim Hoteit, Omar Knio

<figure><img src="image_1.png" /></figure>

## To cite this version:

Mohamad Abed El Rahman Hammoud, Olivier Le Maître, Edriss S Titi, Ibrahim Hoteit, Omar Knio. Continuous and discrete data assimilation with noisy observations for the Rayleigh-Bénard convection: a computational study. Computational Geosciences, 2023, 27, pp.63-79. 10.1007/s10596-022-10180-4. hal-03935869

HAL Id: hal-03935869

https://hal.science/hal-03935869

Submitted on 12 Jan 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Continuous and discrete data assimilation with noisy observations for the Rayleigh-Bénard convection: a computational study

Mohamad Abed El Rahman Hammoud<fnref n="1" /> · Olivier Le Maître<fnref n="2" /> · Edriss S. Titi<fnref n="3" />,<fnref n="4" /> · Ibrahim Hoteit<fnref n="1" /> · Omar Knio<fnref n="5" />

Received: 6 January 2022 / Accepted: 7 November 2022

© The Author(s) 2022

## Abstract

Obtaining accurate high-resolution representations of model outputs is essential to describe the system dynamics. In general, however, only spatially- and temporally-coarse observations of the system states are available. These observations can also be corrupted by noise. Downscaling is a process/scheme in which one uses coarse scale observations to reconstruct the high-resolution solution of the system states. Continuous Data Assimilation (CDA) is a recently introduced downscaling algorithm that constructs an increasingly accurate representation of the system states by continuously nudging the large scales using the coarse observations. We introduce a Discrete Data Assimilation (DDA) algorithm as a downscaling algorithm based on CDA with discrete-in-time nudging. We then investigate the performance of the CDA and DDA algorithms for downscaling noisy observations of the Rayleigh-Bénard convection system in the chaotic regime. In this computational study, a set of noisy observations was generated by perturbing a reference solution with Gaussian noise before downscaling them. The downscaled fields are then assessed using various error- and ensemble-based skill scores. The CDA solution was shown to converge towards the reference solution faster than that of DDA but at the cost of a higher asymptotic error. The numerical results also suggest a quadratic relationship between the $l_2$ error and the noise level for both CDA and DDA. Cubic and quadratic dependences of the DDA and CDA expected errors on the spatial resolution of the observations were obtained, respectively.

**Keywords** Downscaling · Continuous data assimilation · Noisy observations · Rayleigh-Bénard convection

## 1 Introduction

High-resolution representations of the states of dynamical systems are essential in various fields [1]. Downscaling

techniques were introduced to synchronize the fine-scale features of a dynamical system to a reference solution represented by coarse-scale observations [2]. These techniques were applied in various applications, such as oceanic [3, 4] and atmospheric modeling [5, 6]. Downscaling techniques can be categorized into two main approaches: statistical downscaling and dynamical downscaling.

Statistical downscaling establishes a statistical relationship between the observed variables at the coarse-scales and the variables of interest at the fine-scales, based on which, the fine-scales can be estimated [7–10]. Dynamical downscaling methods, on the other hand, constrain the fine-scale model solution to the coarse-scale observations. Nudging techniques are commonly used for dynamical downscaling, where the fine-scale solution is forced towards the coarse-scale data, point-by-point as in grid nudging and low-frequency to low-frequency as in spectral nudging [11]. Several studies compared the performances of statistical and dynamical downscaling methods and reported that both techniques generally perform comparably [6, 12–14].

<figure><img src="image_3.png" /></figure>

<footnote>Omar Knio
omar.knio@kaust.edu.sa</footnote>

<footnote n="1">1 Division of Physical Sciences and Engineering, King Abdullah University of Science and Technology, Thuwal, 23955, Kingdom of Saudi Arabia</footnote>

<footnote n="2">2 Centre de Mathématiques Appliquées, CNRS and Inria, Ecole Polytechnique, 91120, Palaiseau, France</footnote>

<footnote n="3">3 Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge, CB3 0WA, UK</footnote>

<footnote n="4">4 Department of Mathematics, Texas A & M University, College Station, 77843, TX, USA</footnote>

<footnote n="5">5 Division of Computer, Electrical and Mathematical Sciences and Engineering, King Abdullah University of Science and Technology, Thuwal, 23955, Kingdom of Saudi Arabia</footnote>
Continuous Data Assimilation (CDA) is a dynamical downscaling algorithm that constrains the large-scale solution of a system state to coarse-scale observations [15]. The CDA algorithm introduces a nudging term, as a source term in the equations of motion, to compare the coarse scales of the system state to those of the observations. The theoretical grounds of CDA state that the downscaled solution has an exponentially increasing accuracy with time, regardless of the initial conditions [15]. The downscaled solution was also proven to asymptote to the true solution, provided sufficient observations in time. The CDA algorithm was successfully applied to the 2D Rayleigh-Bénard equations [13, 16, 17], 3D Rayleigh-Bénard equations [18], quasi-geostrophic equations [19], 2D Navier Stokes equations [20], 3D Primitive equations [21] and Global Circulation Models (GCMs) [6]. The theoretical framework of CDA with noisy stochastic measurements was more recently established by [22], who showed that the algorithm asymptotes to the true solution with a tolerance limit proportional to the trace of the covariance of the noisy observations.

This work introduces Discrete Data Assimilation (DDA), a downscaling algorithm that is a discrete-in-time counterpart of CDA, arising in spirit of the works of Hayden et al. [23] and Celik et al. [24], where downscaling is performed based on an “implusive” update of the state field, using Fourier-based spectral projections specified in terms of the coarse-scale observations. In contrast, in the present work DDA introduces a source term into the governing equations that is concentrated at the observation times. Thus, whereas CDA continuously nudges the governing equations with the coarse-scale observations of the system state, DDA only nudges the system at the discrete times when observational data are available. To numerically analyze the performances of both downscaling algorithms for the case of noisy coarse-scale observations, we rely on synthetic observational data generated by perturbing a reference solution of the Rayleigh-Bénard convection system in the turbulent regime. These noisy coarse-scale observations of the velocity and temperature are downscaled realization-by-realization using CDA and DDA, which results in a set of downscaled solutions.

The performances of both CDA and DDA were then assessed by analyzing various skill scores. Numerical experiments were also conducted to examine the sensitivity of DDA and CDA to the noise level and temporal and spatial resolutions of the observations. The statistics of the downscaled ensemble of solutions were then examined and the precision of the skill scores estimators were analyzed using a bootstrap analysis. Finally, because velocity observations alone are sufficient to fully describe the downscaled Rayleigh-Bénard solution fields [16], the critical noise level at which noisy temperature observations

should be discarded in favor of employing noisy velocity observations alone was also investigated. The numerical results suggest that the $l_2$-errors resulting from CDA and DDA are proportional to the observational noise variance and have a cubic relation to the spatial resolution. The CDA downscaled fields are also shown to converge to the asymptotic behavior faster than DDA, although it converges to a lower asymptotic error level in comparison to CDA.

The following section provides an overview of the Rayleigh-Bénard convection system, the CDA and DDA algorithms, and the skill scores to assess their performances. Section 3 describes the experimental setting, including the generation of noisy observational data, and the model and downscaling parameters. The results of the numerical experiments are then presented in Section 4. The main conclusions of the study are outlined in Section 5.

# 2 Preliminaries

## 2.1 Rayleigh-Bénard convection system

The Rayleigh-Bénard convection is a thermo-fluidic instability driven by the temperature difference between a hotter bottom boundary and a cooler top boundary. At high Ra, chaos dominates the flow with no particular structures to be observed [25, 26]. The system under study consists of a 2D rectangular periodic channel of vertical height $L_y$ and horizontal period $L_x$. The 2D Boussinesq equations are then recast into their non-dimensional form, using the same scaling used in [27, 28], such that the equations are solved in a fundamental spatial domain $\Omega = [0, L_x/L_y] \times [0, 1]$. The solution of the governing equations are the dimensionless velocity vector $\mathbf{u} = (u, v)$, pressure $p$ and temperature $\Theta$. These equations are expressed as:

$$ \nabla \cdot \mathbf{u} = 0, \qquad (1) $$

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} + \nabla p = \frac{Pr}{\sqrt{Ra}} \nabla^2 \mathbf{u} + Pr \Theta \mathbf{e}_y, \qquad (2) $$

$$ \frac{\partial \Theta}{\partial t} + (\mathbf{u} \cdot \nabla) \Theta = \frac{1}{\sqrt{Ra}} \nabla^2 \Theta + \mathbf{u} \cdot \mathbf{e}_y, \qquad (3) $$

where $Ra := (g\alpha\Delta T L_y^3)(\nu\kappa)^{-1}$ is the Rayleigh number, $g$ is the gravitational acceleration, $\alpha$ is the thermal expansion coefficient of the fluid, $\Delta T$ is the dimensional temperature difference between the top and bottom boundaries, $\nu$ is the kinematic viscosity, $\kappa$ is the thermal diffusivity coefficient, $Pr = \nu\kappa^{-1}$ is the Prandtl number, $t$ is time and $\mathbf{e}_y$ is the unit basis vector in vertical direction.

This system of equations can be solved once the initial and boundary conditions are specified. The temperature and velocities are initialized with independent, spatially
uncorrelated random noise sampled point-by-point from uniform distributions $U(-0.2, 0.2)$ and $U(-0.1, 0.1)$, respectively, where $U(a, b)$ refers to a uniform distribution on the interval $(a, b)$. Periodic boundary conditions were specified in $x$, whereas at the horizontal boundaries, no-slip and no-penetration along with isothermal conditions were enforced. Specifically, the boundary temperatures are set as $\Theta(t; x, 0) = \Theta(t; x, 1) = 0$.

The system of equations were solved numerically using finite differences on a staggered grid using the conservative form of the governing equations. A uniformly spaced staggered grid with $n_x \times n_y$ cells and spacing $\delta x = \delta y = h$ was adopted. The size of the computational cell is $\delta x = L_x/(L_y \times n_x)$ horizontally and $\delta y = 1/n_y$ vertically, where $n_x$ and $n_y$ are the number of horizontal and vertical grid points, respectively. A third-order Adam-Bashforth scheme was employed using a fixed time step $\delta t$. A second-order central differencing scheme was selected for the advection, diffusion and source terms. Pressure-projection based on the fractional step method was utilized to satisfy the incompressibility condition, which is manifested in the conservation of mass equation. Fast Fourier Transform (FFT) based algorithms were implemented to efficiently solve the parabolic and elliptic systems.

## 2.2 Continuous and discrete data assimilation

CDA is a nudging-based downscaling algorithm that employs coarse-scale observational data to constrain the large-scale features of the model solution in order to dynamically recover the fine-scale features. This is achieved by introducing a source term, proportional to the discrepancy between the coarse scales of the model solution and those of the observational data, in the equations of motion and energy [15]. CDA is applied here to the 2D Rayleigh-Bénard convection for downscaling coarse observational data by solving the following system of equations:

$$ \nabla \cdot \mathbf{w} = 0, \qquad (4) $$

$$ \frac{\partial \mathbf{w}}{\partial t} + (\mathbf{w} \cdot \nabla) \mathbf{w} + \nabla p = \frac{Pr}{\sqrt{Ra}} \nabla^2 \mathbf{w} + Pr \Psi \mathbf{e}_y + \mu_u (I_{h^o} (\mathbf{u}^o) - I_{h^o} (\mathbf{w})), \quad (5) $$

$$ \frac{\partial \Psi}{\partial t} + (\mathbf{w} \cdot \nabla) \Psi = \frac{1}{\sqrt{Ra}} \nabla^2 \Psi + \mathbf{w} \cdot \mathbf{e}_y + \mu_{\Theta} (I_{h^o} (\Theta^o) - I_{h^o} (\Psi)), \quad (6) $$

where $\mathbf{w} = (\tilde{u}, \tilde{v})$ is the downscaled velocity vector, $\Psi$ is the downscaled temperature, $\mu = (\mu_u, \mu_\Psi)$ are non-negative constants called the nudging parameter, and $\mathbf{u}^o$ and $\Theta^o$ are the coarse-scale observations of the velocity and temperature, respectively. Note that the downscaling system is solved using the same boundary conditions as the reference system. Finally, $I_{h^o}$ is an interpolation

operator of an interpolant function $\phi$ for uniformly spaced measurements at a distance $h^o$ in $\Omega$, which is expressed as follows:

$$ I_{h^o}(\phi(\mathbf{x})) = \sum_{k=1}^{N_{h^o}} \phi(\mathbf{x}_k) \chi_{Q_k}(\mathbf{x}), \quad (7) $$

where $\mathbf{x} = (x, y)$, $Q_k$ are disjoint subsets such that $\text{diam}(Q_k) \le h^o$ for $k = 1, \dots, N_{h^o}$, $N_{h^o}$ is the number of observational points, $\bigcup_{j=1}^{N_{h^o}} Q_j = \Omega$, $x_k \in Q_k$ and $\chi_E$ is the characteristic function of the set $E$.

It was proven that regardless of the initial conditions used for $\Psi$ and $\mathbf{w}$, the CDA solution converges exponentially to the true solution under the conditions [15] outlined. Specifically, the system converges if sufficiently enough coarse observational data that is spaced at a sufficiently small $h^o$ are provided and the $\mu$'s are appropriately selected. CDA, thus, ensures that, for dissipative systems, the downscaled solution is synchronized with the reference solution for all times beyond the last observation time.

Inspired by the works of [15, 23] and [24], we propose the DDA algorithm in which nudging is performed only at discrete instances when observations are available, while solving the original system at all other steps. Specifically, DDA simulates (5) and (6) when observations are available, and integrates the original (2) and (3) at all other times. The discrete forms of the energy equations for CDA and DDA are presented in order to provide additional insight on the implementation of the proposed downscaling algorithms as:

$$ \frac{\Psi^{n+1} - \Psi^n}{\Delta t} + (\mathbf{w} \cdot \nabla \Psi)^n = \frac{1}{\sqrt{Ra}} \nabla^2 \Psi^n + (\mathbf{w} \cdot \mathbf{e}_y)^n + \mu_{\Theta} [I_{h^o}(\Theta^o(t_i^o)) - I_{h^o}(\Psi^n)] \quad (8) $$

for CDA and

$$ \frac{\Psi^{n+1} - \Psi^n}{\Delta t} + (\mathbf{w} \cdot \nabla \Psi)^n = \frac{1}{\sqrt{Ra}} \nabla^2 \Psi^n + (\mathbf{w} \cdot \mathbf{e}_y)^n + \sum_{i=1}^{n_{y^o}} \delta(t - t_i^o) \mu_{\Theta} [I_{h^o}(\Theta^o(t_i^o)) - I_{h^o}(\Psi^n)] \quad (9) $$

for DDA, where $\delta$ is the dirac delta function. In (8) and (9) the algorithm advances from time step $n$ to $n+1$ and $t_i^o$ is the instance of the $i^{th}$ observation. CDA, therefore, applies the nudging term that employs the $i^{th}$ observation between $t_i^o$ and $t_{i+1}^o$. On the other hand, DDA includes the nudging term only at the time step when the observation was made and solves the original system of equations otherwise.

The assimilation of perfect observations was analyzed both theoretically and numerically for CDA. The effect of noisy observations on the convergence rate and the asymptotic error level have only been examined theoretically by [22]. Bessah et al. [22] showed that the expected $\ell_2$-norm of the error, $\Lambda$, between the downscaled ($Z_d$) and true ($Z_t$)
solution fields is proportional to the trace of the covariance of the observational noise, which is expressed as follows:

$$ \Lambda := \mathbb{E} \left[ \int (Z_d - Z_t)^2 d\Omega \right] \propto C \cdot \mathrm{Tr} [\mathrm{Cov}] \quad (10) $$

where Z is an arbitrary variable, dΩ represents a differential area in the domain, C is a functional that depends on the numerical and downscaling parameters, and Cov is the covariance of the noise. One objective of this work is to numerically investigate the asymptotic error level and the convergence rate. Another objective is to estimate the functional dependence of Λ on the covariance of the noise, and spatial and temporal resolutions of the observations, based on which the CDA and DDA algorithms were compared.

## 2.3 Generation of random observations

The solution of (1)-(3) with the initial and boundary conditions described in Section 2.1 is considered as the reference solution. Coarse observations of the system states were generated by subsampling points every h° spatially and δt° temporally from the reference solution. (Note that as further discussed below, the reference solution is obtained on a sufficiently fine computational grid that it can be practically treated as the exact solution.) The spatial and temporal resolutions of the observational grid are characterized using the non-dimensional ratios, called spatial and temporal downscaling factors, defined as $R = h°/h$ and $S = \delta t°/\delta t$, respectively.

The coarse-scale observations of the reference solution are then perturbed by adding randomly sampled independent and identically distributed (iid) white noise sampled from a Gaussian distribution with zero mean and standard deviation σ. In this work, σ = (σT, σu) is also referred to as the noise level, where the noise of the two variables are independent to mimic different measurement sources, a practical scenario encountered with observational data. Each noisy observed realization is then downscaled separately to generate an ensemble of downscaled solutions, each representing a possible solution trajectory of the system. The performance of the downscaling algorithms was then evaluated quantitatively based on the asymptotic error level between the downscaled and reference solutions, and the rate at which they reach the asymptotic behavior, i.e. the convergence rate.

## 3 Experimental setup

The experiments were conducted using the Rayleigh-Bénard convection with a high Rayleigh number of Ra = 2 × 10⁸ and Pr = 0.7, which corresponds to a highly

chaotic system that mimics environmental flows [29, 30]. Numerically, the system of equations are solved on a uniformly spaced staggered grid composed of 1200 × 400 grid points (for scalar fields) discretizing a domain of length Lx = 3 and height Ly = 1. The time step δt was set to 5 × 10⁻⁴, leading to a maximum Courant number of about 0.15 and a maximum grid Reynold number of O(10). This ensures that fine scale dynamics are suitably captured. In other words, the computed solution is essentially insensitive to further refinement of the spatial or temporal resolution of the computational grid.

The sensitivity of the convergence rate and the asymptotic error level, of both CDA and DDA, to the downscaling parameters and the observation noise level were investigated to determine the corresponding functional dependencies. Specifically, the impact of the observation noise level was numerically studied by analyzing the sensitivity of the downscaled solution statistics to different values of σ. The solution variables were initialized with a random field and the computations were carried out over a sufficiently long interval to observe convergence. This generally occurred within t = 25 with DDA, and for a much shorter window (t = 3) with CDA. Moreover, the asymptotic error level was analyzed for the final steps, which, in this study, correspond to $t_f^{DDA} = 49.9$ for the DDA and $t_f^{CDA} = 15$ for the CDA.

A 0th order interpolant operator was selected because it is computationally cheapest with no noticeable impact on the convergence rate, as was also suggested by [13]. The nudging factor for both the energy and momentum equations were chosen the same with μ° = μu = μ to limit the number of hyper-parameters. In the present work, a grid search was performed to find the best μ, where the downscaled fields were ensured to converge to the reference solution. Note that the theory in [15] predicts the existence of a suitable range of parameters, and computational experience [13] suggests that CDA predictions always converge so long as μ is selected within this range. The performance of the downscaling algorithms also depends on the spatial and temporal resolutions of the available observations. The sensitivity of both CDA and DDA to these parameters were, therefore, also investigated.

The quality of the downscaled solutions were assessed by computing various skill scores to quantify the errors relative to the reference solution and variability relative to the mean solution. The skill scores were computed for each of the solution variables and included the root-mean-squared-error, which is expressed as:

$$ RMSE(t) = \sqrt{\frac{1}{n_x n_y} \sum_{i=1}^{n_x n_y} [Z_d(x_i; t) - Z_t(x_i; t)]^2}, \quad (11) $$
where $x_i = (x_i, y_i)$ is the position of the grid point and $t$ is time; the relative root-mean-squared error:

$$RRMSE(t) = \frac{\sqrt{\sum_{i=1}^{n_x n_y} [Z_d(x_i; t) - Z_t(x_i; t)]^2}}{\sqrt{\sum_{i=1}^{n_x n_y} [Z_{true}(x_i; t)]^2}}. \quad (12)$$

and the absolute error:

$$AE(t) = \frac{1}{n_x n_y} \sum_{i=1}^{n_x n_y} |Z_d(x_i; t) - Z_t(x_i; t)|. \quad (13)$$

The variability of the downscaled solutions was evaluated as the average ensemble spread as follows:

$$AES(t) = \sqrt{\frac{1}{N_e - 1} \sum_{k=1}^{N_e} \sum_{i=1}^{n_x n_y} (Z_d^k(x_i; t) - \bar{Z}_d(x_i; t))^2}, \quad (14)$$

where $N_e$ is the number realizations forming the ensemble of downscaled solutions and $\bar{Z}_d(x_i; t)$ is the average of the samples at $x_i$ and time $t$. Finally, the measure introduced by [22] on the expected $\ell_2$-error of the downscaled members, denoted by $\Lambda$, was also computed based on the definition in (10) for the last time step of the simulation. This enables us to numerically assess the theoretical results and to further investigate the functional dependence on downscaling parameters.

**Remark** It is instructive to contrast the resolution of the coarse observation grids to those of the integral time scales of the flow. Based on the definitions above, the spatial and temporal resolutions of the observations grid are given by $h^0 = \mathcal{R}h$ and $\delta t^0 = S\delta t$. In the downscaling experiments conducted below, the downscaling ratios are selected such that $5 \le \mathcal{R} \le 25$ whereas $1 \le S \le 50$. With these choices, $h^0$ and $\delta t^0$ remain substantially smaller than the flow integral space scale and time scale, respectively. In particular, the dimensionless large-scale turnover time, $\tau \equiv 2L_y/u_{max}$, estimated based on the reference solution, is approximately 3, which indicates that $\delta t^0/\tau$ remains much smaller than 1 for the entire range of the $S$ values considered.

# 4 Results

This section presents and discusses the results of the numerical experiments illustrating the implementation of the DDA and CDA algorithms to downscale coarse scale noisy observations of the Rayleigh-Bénard flows. The DDA results are mostly presented in the main text and when not included the CDA results are mirrored in the supplementary material, but are analyzed in the text for direct comparison.

## 4.1 Sensitivity to noise levels

The sensitivity of each of the CDA and DDA algorithms to observation noise levels were examined by first selecting a nominal case, defined by a nominal noise level equal to 10% of the range of the solution variable at the final time step. The noise levels for the nominal case were set to $(\sigma_T, \sigma_u) = (0.1, 0.05)$. Higher and lower noise levels are then considered by scaling both noise levels simultaneously. Consequently, the results presented indicate $\sigma_T$ only.

Noisy realizations of the coarse-scale observations of the system are sampled as described in Section 2.3 and are then downscaled using CDA and DDA, respectively. The DDA and CDA nudging parameters are set to 7.0 and 3.0, respectively, for an observational grid resolution of $S = 10$ and $\mathcal{R} = 5$. The nudging factor was tuned by solving the downscaling equations for multiple values of $\mu$, and the solution with lowest RRMSE at the final time step is used throughout. Downscaling the noisy realizations was conducted starting from a random field to verify that, given sufficient observations and appropriate nudging, both algorithms converge regardless of the initial conditions.

Figure 1 presents the time evolution of the AE, RMSE, RRMSE of the DDA downscaled temperature solutions. The figure also illustrates the AES for ensembles composed of 10, 20, 25, 30, 40 and 50 downscaled solutions. The evolution of the RRMSE from the CDA simulation is also shown for comparison. The results of the velocity components are similar to those of the temperature; they are briefly discussed in this section with the corresponding results provided in the Supplementary.

The time evolution of both the AE and RMSE of the DDA downscaled solutions for all 50 members of the ensemble are illustrated in Fig. 1a and b, respectively. The skill scores decrease in time until reaching a plateau at approximately $t = 30$. Figure 1c shows a rapid decay of the AES, for all considered ensemble sizes, following similar patterns to those of AE and RMSE. The plots also indicate that the smallest AES values are at the plateau. Similar to the downscaled temperature fields, the errors of the downscaled velocity fields also decay to an asymptotic error level, and achieve AES values that are comparable to those of the downscaled temperature fields.

Figure 1d and e illustrate the decay of the RRMSE for both DDA and CDA downscaled fields. These plots show that with CDA the RRMSE decays smoothly and exponentially in time, whereas with DDA the RRMSE decays non-monotonically and at a slower rate. As time evolves, the skill scores of both algorithms asymptote to a plateau, with the asymptotic error skill scores corresponding to DDA being lower than that of CDA. Specifically, the asymptotic RRMSE of the DDA downscaled fields were approximately three times smaller than that corresponding
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 1** Evolution of different skill scores corresponding to the DDA downscaled temperature fields, each member is driven by observations of the nominal noise level of ($\sigma_T, \sigma_u$) = (0.1, 0.05). Each plot contains curves corresponding to 50 downscaled samples, the AES is presented for different number of realizations as indicated</figcaption></figure>

to CDA, and is approximately one third of the observations' noise level. The RRMSE of the downscaled velocity components experience a similar evolution to that of the downscaled temperature fields, where the RRMSE decays smoothly and exponentially with CDA and non-monotonically with DDA. This error behavior can be attributed to the nature of the nudging involved with each algorithm. Specifically, CDA continuously nudges the system of equations, which introduces observational noise at every integration time step. On the other hand, DDA nudges the system of equations at discrete times, allowing the unforced system of equations to relax this sudden change imposed by the nudging at discrete time steps. Hereafter, we only discuss the evolution of the RRMSE of the downscaled temperature solution because the skill scores evolve similarly in time.

The reference temperature solution along with its difference with the DDA downscaled solutions are presented at different time instants in Fig. 2. The snapshots first show the initial condition of the reference system and the random initialization of the DDA systems, which is expressed by the noisy error field. The solution at the intermediate time shows a decrease in the errors across the domain, mostly concentrated near the plumes, and are of the order of the observation noise level. At a later time, the errors practically vanish across the domain, with the highest error values being an order of magnitude smaller than the noise level. The highest errors occur near detachment and mixing zones due to the more complex fine-scale dynamics, as expected.

To provide insight into the effect of the noise level on the downscaled fields, the reference temperature field and the difference of the DDA downscaled temperature fields with this reference are presented for different snapshots in Fig. 3. In all cases, the results indicate that the difference between the downscaled and reference temperature fields is largest around plumes and their corresponding boundaries. At an intermediate time ($t = 15$), the discrepancy fields reveal structures whose shapes and amplitudes depend on the observation noise level, especially when the latter is large. This is not surprising because in this situation the downscaled solution is still far from the reference solution. As further observations are assimilated in time ($t = 45$), the impact of increasing the noise level essentially manifests itself as an amplitude scaling of the discrepancy field. Thus, when the error plateau is reached, the discrepancy fields exhibit similar structure, with amplitude differences that scale with the observation noise level.

The distributions of the asymptotic skill scores corresponding to the DDA downscaled temperature field at the final time step are shown in Fig. 4. The box plots in Fig. 4a and b depict an increase in the asymptotic AE and RMSE quadratically with the noise level. Figure 4c shows the decrease in the AES for increasing ensemble size at a rate of
<figure><img src="image_9.png" /><figcaption>(a) t = 0</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) t = 10</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) t = 45</figcaption></figure>

<caption>**Fig. 2** Instantaneous snapshots of the temperature distribution for the noise-free Rayleigh-Bénard simulation along with the difference between the downscaled sample solutions and the truth. Snapshots are</caption>

shown for the initial time step (top), an intermediate time step (middle) and a time step when the algorithm had converged (bottom)
<figure><img src="image_12.png" /><figcaption>**Fig. 3** Instantaneous snapshots of the temperature distribution for the noise-free Rayleigh-Bénard simulation along with the difference between the downscaled solutions and the truth. Snapshots are shown for an intermediate time step (top) and a time step when the algorithm had converged (bottom). At both instants, plots are generated for different noise levels, as indicated</figcaption></figure>

(a) $t = 15$

<figure><img src="image_13.png" /></figure>

(b) $t = 45$

approximately $1/\sqrt{N_{ens}}$, as expected for random sampling of the observation noise. The figure also illustrates a linear decrease of the AES with decreasing noise level and a fixed ensemble size. Finally, the curves in Fig. 4d present the increase of $\Lambda$ with $\sigma$ on a log plot for both CDA and DDA. The results suggest that $\Lambda$ is proportional to the variance of the noise for both the CDA and the DDA downscaled fields, consistent with the theory of [22].

## 4.2 Sensitivity to the observation's temporal sampling rate

This section analyzes the impact of the temporal observation sampling rate on the DDA and CDA solutions. Numerical experiments were performed with different observation frequencies by systematically varying $S$ within the broad range $1 \le S \le 50$. In doing so, the noise level was held fixed with $\sigma = (0.1, 0.05)$ and $R = 5$.

Figure 5 shows the time evolution of the temperature RRMSE resulting from DDA for $S = 1, 10$ and $24$. For $S = 1$, the RRMSEs decay monotonically to a plateau. For $S = 10$ and $24$, however, the behavior becomes non-monotonic. In both cases, the RRMSE initially drops

rapidly, then increases before decreasing again to eventually reach a plateau. As $S$ increases, the downscaled fields exhibit a wider spread. Specifically, for $S = 1$, the RRMSE curves are almost indistinguishable, whereas for larger $S$, the variability in the downscaled fields becomes substantial, though a smaller mean error is achieved for $S = 10$. CDA also seems to be less sensitive to $S$ compared to DDA, where the time evolution of the temperature RRMSE of different members is almost identical and the plateau is reached almost immediately.

Figure 6 shows the distribution of the asymptotic values of AE, RRMSE and AES for the DDA downscaled temperature fields with varying $S$ at the final time step. The behavior of $\Lambda$ in response to varying $S$ is presented for both CDA and DDA for comparison. Specifically, Fig. 6a and b show the box plots describing the AE and RRMSE of the temperature corresponding to each of the considered observation frequencies. For $S \le 15$, the distributions of the AE and RRMSE are tight with little variability. The variability of these skill scores increases as $S$ increases. The box plots also show that both the AE and RRMSE initially decrease with increasing $S$ until reaching their minimum at $S = 10$, after which, the errors increase nonlinearly
<figure><img src="image_14.png" /><figcaption>(a) Absolute Error</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) RRMSE</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) AES</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(d) Λ</figcaption></figure>

**Fig. 4** Distribution of the skill scores (indicated) based on the 50 sample solutions for the ensemble of temperature solutions of the (a)-(d) DDA solution, and (d) CDA solution at their respective final time steps. Sample solutions were generated by varying the noise levels (indicated by the horizontal axis), while keeping the remaining parameters fixed

<figure><img src="image_18.png" /><figcaption>(a) $S = 1$</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) $S = 10$</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(c) $S = 24$</figcaption></figure>

**Fig. 5** Evolution of the RRMSE of 50 individually downscaled temperature solutions for: (a) $S = 1$, (b) $S = 10$ and (c) $S = 24$. Note that with $S = 24$ the simulations were carried out till $t = 120$, as reflected in the inset

with increasing $S$. In particular, for $S \ge 25$, both AE and RRMSE show large values indicating a substantial deviation of the downscaled solution from the reference solution. Figure 6c plots the variation of AES for different $S$ and number of downscaled fields. As expected, the results suggest that for a fixed $S$, the AES decreases as $1/\sqrt{N_{ens}}$. For a fixed given number of downscaled fields, the AES decreases to a minimum as $S$ increases from 1 up to 10, beyond which the AES continues to increase following a similar behavior to that of the AE and RRMSE. Finally, Fig. 6d illustrates the variation of $\Lambda$ with $S$ for both the CDA and DDA downscaled temperature fields. $\Lambda$ varies linearly with $S$ for CDA and exhibits larger values than those obtained with DDA, except for extreme values of $S$.
<figure><img src="image_21.png" /><figcaption>(a) Absolute Error</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(b) RRMSE</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(c) AES</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(d) Λ</figcaption></figure>

**Fig. 6** Frames (a) and (b) show the box plots of the temperature skill scores, AE and RRMSE respectively, for 50 downscaled solutions using DDA, at the final time step. Frame (c) depicts the AES for different number of realizations of the DDA downscaled fields. Frame (d) illustrates the dependence of $Λ$ on $R$ for DDA and CDA as indicated by the legend

On the other hand, $Λ$ evolves non-monotonically with $S$ for the DDA downscaled fields, and decreases to a minimum as $S$ increases from 1 up to 10 before nonlinearly increasing with increasing $S$.

The behavior of DDA errors is attributed to the discrete nature of the nudging, which allows the solution to evolve based on the exact system of equations between the nudging time steps. At very small values of $S$, however, nudging is frequently applied, and thus, the solution becomes increasingly affected by the observation noise leading to an additional error.

## 4.3 Sensitivity to the observation's spatial resolution

The CDA was theoretically proven to converge if the observations are not too far apart [15]. Specifically, the resolution of the observational grid must fall within the range of the dissipation spectrum. In this experiment, the response of CDA and DDA to different spatial resolutions of the observational grid is examined. Downscaling runs were then performed for $R$ of 5, 15 and 25, all with $S = 10$. The joint effect of coarse observations in space and time is also examined with $S = 25$ and $R = 5$ and 15. The nudging factors of DDA and CDA were set to 7.0 and 3.0, respectively. The observational noise levels are set to their nominal levels, $(σ_T, σ_u) = (0.1, 0.05)$.

Figure 7 shows the time evolution of the temperature RRMSE for the different combinations of $S$ and $R$ values. In all cases, the RRMSE exhibits an initial rapid drop followed by an increase before slowly settling down to a plateau. For $S = 10$, the height of the plateau noticeably increases as $R$ increases. Specifically, for $R = 15$, all the RRMSEs of the downscaled realizations experience an initial rapid drop followed by a small increase ending with a plateau. As $R$ increases, however, the RRMSE curves clearly exhibited a larger spread and reached a higher value at the plateau, which took a longer period to reach. On the other hand, for $S = 25$ and for both $R$ values examined, the RRMSE curves exhibit a smaller rise in the plateau with increasing $R$ in comparison to the case of $S = 10$. The monotonic decrease in the CDA's RRMSE is lost with increasing $R$, as seen in the corresponding results in the Supplementary material. In particular, for $R = 15$ and 25, the RRMSE rapidly decreases, then increases for a short period before decreasing to a final plateau. The DDA seems to be more sensitive to $R$ for the smaller value of $S$ indicated by the increase in the value of the plateau. These results suggest that the coarser the spatial resolution of the observational grid is, the less accurate the downscaled solution becomes. A thorough assessment of the behavior of the skill scores at the final time step is presented in the following.

<figure><img src="image_25.png" /></figure>
<figure><img src="image_26.png" /><figcaption>**Fig. 7** Evolution of the RRMSE of 50 downscaled temperature fields corresponding to DDA. Plots are generated for different $\mathcal{R}$ and $\mathcal{S}$ as indicated</figcaption></figure>

Figure 8 presents the box plots of AE, RRMSE and AES for $\mathcal{S} = 10$ and the different tested $\mathcal{R}$ values for the DDA downscaled temperature fields at the last time step. The figure also displays the behavior of $\Lambda$ for both CDA and DDA downscaled temperature fields at the last time step. The plots suggest a cubic increase of the means of both the AE and RRMSE and an increase in the solutions' spread with increasing $\mathcal{R}$. The AES increases linearly with increasing $\mathcal{R}$ and decreases with the square-root of the number of downscaled realizations. Finally, Fig. 8d shows the variation of $\Lambda$ as a function of $\mathcal{R}$ for both the CDA and DDA downscaled fields. Specifically, $\Lambda$ increases with increasing $\mathcal{R}$ for both algorithms. $\Lambda$ increases quadratically with $\mathcal{R}$ with CDA, whereas with DDA, $\Lambda$ is a cubic function of $\mathcal{R}$. For comparison, Fig. 9 shows the box plots of the AE and RRMSE at the final time step for the case with $\mathcal{S} = 25$ and the considered values of $\mathcal{R}$. Both AE and RRMSE increase when $\mathcal{S}$ is larger and increase for increasing $\mathcal{R}$. In comparison to the results of Fig. 8, the mean of the asymptotic errors almost triple, for a fixed $\mathcal{R}$, when $\mathcal{S}$ is larger. In addition, the variance of the downscaled fields becomes larger with larger $\mathcal{S}$.

## 4.4 Distribution of the downscaled fields

Here, we analyze the distribution of the downscaled fields by applying a statistical test. In this experiment, the noise level was set to the nominal value, the DDA and CDA nudging parameters were fixed to 7.0 and 3.0, respectively, with $\mathcal{S}$ and $\mathcal{R}$ equal to 11 and 5, respectively. 200 solution samples were then generated and their profiles in the planes $y = 0.25, 0.50$ and $0.75$ were selected. The Kolmogorov-Smirnov test was then applied to the solution points in these planes at discrete times to assess their distribution. The goal is to test for Gaussianity of the downscaled solutions in space at different instances in time.

Figure 10 plots snapshots of the p-value and hypothesis test for $t = 30$ and $t = 49.9$, representing an intermediate and the final time steps. The plots for $t = 30$ show several locations where the hypothesis tests positive, hypothesis is equal to 1 or equivalently p-value less than 0.05, meaning that the samples fail the Gaussianity test. However, the plots at $t = 49.9$ indicate that only two out of the total 3600 tested locations fail the Gaussianity test, suggesting that, for iid Gaussian noise, the downscaled solutions are normally
<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /><figcaption>**Fig. 8** Frames (a) and (b) show the box plots of the temperature skill scores, AE and RRMSE respectively, for 50 downscaled solutions using DDA with $S = 10$ and the indicated $\mathcal{R}$, at the last time step. Frame (c) depicts the AES for different number of realizations. Frame (d) illustrates the dependence of $\Lambda$ on $\mathcal{R}$ for both DDA and CDA</figcaption></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /><figcaption>**Fig. 9** Box plots of the temperature skill scores: (a) AE and (b) RRMSE for 50 downscaled solutions using DDA with $S = 25$ and the indicated $\mathcal{R}$</figcaption></figure>

distributed once the plateau is reached. Similar results were obtained for the CDA downscaled samples, but those passed the Gaussianity test for all locations at an earlier time due to its faster convergence rate.

## 4.5 Precision of the skill scores and a lower error solution

In these experiments, we first assess the precision of the estimate of the RRMSE by means of a bootstrap methodology. Second, we explore the potential of using the average of the downscaled solutions as the best solution. To assess the precision of the skill score estimates, bootstrap sampling is applied to subsets of various sizes of the ensemble of 200 downscaled solutions of the previous section. A total of 500 bootstrap subsets were generated for each skill score.

The average of each bootstrapped skill score subset is considered as a bootstrap estimate of the skill score. These samples were then input to a Kernel Density Estimator (KDE) to generate the probability density function (PDF) of the skill scores.

Figure 11 presents the resulting distributions for subsets of different sizes. Specifically, Fig. 11a illustrates the distribution of the average RRMSE estimate for the
<figure><img src="image_33.png" /><figcaption>(a) Hypothesis Test, $t = 30$</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>(b) Hypothesis Test, $t = 49.9$</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>(c) P-Value, $t = 30$</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(d) P-Value, $t = 49.9$</figcaption></figure>

**Fig. 10** Plots of the spatial distribution of the hypothesis value and p-value for the Kolmogorov-Smirnov test. In each plot, distributions are generated based on horizontal profiles at $y = 0.25$, $y = 0.5$ and $y = 0.75$ as indicated. Plots are generated at $t = 30$ and $t = 49.9$ using a 200 realizations of the downscaled temperature solutions

different subsets and suggests that the variance in the estimated average RRMSE decreases with subset size. The box plots describing the distribution of the boostrapped

RRMSE are illustrated in Fig. 11b. The figure shows a shrinking variability of the estimated average RRMSE with increasing subset size. In particular, by bootstrapping a subset with as few as 20 realizations, the standard deviation of the distribution of the skill scores estimate can be approximately halved in comparison to the subset of 10 realizations. Moreover, beyond 30 solution samples, the standard deviation of the average RRMSE asymptotes to a steady value.

Because the observational noise is unbiased, one may anticipate the mean solution to provide a better estimate of the reference solution than the individual realizations. The skill scores associated with the ensemble mean were computed for different subset sizes and compared to the distribution of the skill scores obtained from the 200 nominal case solutions. Figure 12 illustrates the box plots of the asymptotic temperature RRMSE for both CDA and DDA downscaled fields, along with the curve of the RRMSE of the average solution for different numbers of realizations. The results show that, for both the CDA and DDA and any number of realizations, the RRMSE of the average solution is lower than the RRMSE of any downscaled realization. The RRMSE of the average solution is also seen to decrease with larger number of realizations. Specifically, the RRMSE of the DDA downscaled solution dropped from approximately 3.5% to below 1% for more than 20 realizations. Moreover, comparing the distributions from CDA to those from DDA, the results suggest that the DDA yields errors three times smaller than those of the CDA downscaled solutions. The CDA, however, shows a tighter distribution with a smaller spread between the downscaled fields.

## 4.6 Observing a noisy solution

In a practical setting, only a single noisy observation of the state is usually available. Here, we generate samples of noisy observations starting by downscaling the given noisy coarse-scale observations to obtain a high-resolution downscaled reference. Samples of noisy observations were then generated by adding iid noise sampled from a Gaussian distribution to this downscaled reference solution. The resulting noisy observations are then considered as perturbed noisy observations. The sampled noise, at the nominal level, is added to the downscaled reference solution, as in Section 4.1. By doing so, we investigate how CDA and DDA behave with the sampled noisy observations compared to the case when realizations of noisy observations of the truth are readily available. Both CDA and DDA algorithms were then used to downscale these perturbed noisy observations for $S = 10$ and $R = 15$. The downscaled fields are compared to the reference fields to assess the performance of the downscaling algorithms.
<figure><img src="image_37.png" /><figcaption>**Fig. 11** Distribution of the RRMSE of the bootstrapped temperature solutions for different number of realizations, as indicated, at the last time step. Plot (a) shows the PDF of the bootstrapped RRMSE of the temperature, generated using KDE, for different number of realizations as indicated, and (b) the box plots representing the distribution of the 500 bootstrap RRMSE samples for different number of realizations as indicated</figcaption></figure>

Figure 13 presents the box plots describing the distribu-
tion of the RRMSE of the DDA downscaled temperature
fields relative to the reference temperature field at the final
step. The box plots illustrate that the RRMSE increases by
20% when the perturbed noisy observations were assimi-
lated. Likewise, the variance of the downscaled solutions
also increases and the RRMSE of the average solution is
amplified by approximately 3.5 times when downscaling
was performed on the perturbed noisy observation. Never-
theless, the asymptotic RRMSE of the average solution was
smaller than any of the asymptotic RRMSEs of the down-
scaled solutions, suggesting that it is always worth down-
scaling the perturbed noisy observations if computational
resources are available.

4.7 Relevance of temperature observations

The CDA was shown to converge to noise-free observational
data using coarse-scale velocity observations alone [16].
The convergence of the algorithm, however, is faster when
both temperature and velocity observations are assimilated
[13]. Here, we assess the threshold for the temperature

<figure><img src="image_38.png" /><figcaption>**Fig. 12** (a) Boxplots describing the distribution of the RRMSE of the downscaled temperature solutions for the DDA solution at the last time step (blue) and CDA solution at *t* = 15 (orange) for different sample sizes. (b) RRMSE curves of the average solution, at the final time, of the DDA (blue) and the CDA (orange) downscaled fields for different numbers of downscaled samples</figcaption></figure>

noise level beyond which temperature observations become
disruptive to the convergence of CDA and DDA, and thus
should be discarded.

In this experiment, observations of the velocity and
temperature were provided for $S = 15$ and $R = 11$
with a fixed velocity noise level of $\sigma_u = 0.05$ and three
temperature noise levels: $\sigma_T = 0.05, 0.1$ and $0.15$. An
experiment was also conducted by downscaling using the

<figure><img src="image_39.png" /><figcaption>**Fig. 13** Box plots comparing the distributions of the RRMSE of the temperature sample solutions for the cases where observational data is provided based on the (left) truth and (right) noisy downscaled fields. The blue crosses represent the RRMSE of the average of all the downscaled samples</figcaption></figure>

<figure><img src="image_40.png" /></figure>
<figure><img src="image_41.png" /><figcaption>Fig. 14 Boxplot of the RRMSE of the DDA downscaled temperature solutions for temperature observations with low, medium, high temperature noise level, and no temperature observation. Each box plot is generated from a 50 downscaled solutions at the last time step</figcaption></figure>

coarse velocity observations only, which is equivalent to setting $\mu_{\ominus} = 0$.

The time evolution of the RRMSE of the downscaled temperature fields illustrates that, in all tested cases, the RRMSE decreases until it reaches a plateau at approximately $t = 40$. All RRMSE curves suggest similar characteristics before reaching the plateau, however, the spread is larger in the case when temperature was not assimilated. Figure 14 presents the box plots describing the distribution of the asymptotic RRMSE of the DDA downscaled temperature. The downscaled solution without assimilating temperature observations achieves lower RRMSE than the case of high temperature noise level, comparable RRMSE with the case of medium temperature noise, and a larger RRMSE than that of the small temperature noise level. The case with no temperature observations shows the largest variance compared to all tested cases. This suggests that, in the case of medium to high noise, it might be better to discard temperature observations and opt for solely assimilating medium-noise velocity observations. On the other hand, if reliable temperature observations are available, it would be useful to jointly downscale with those of the velocity to obtain a more accurate solution.

## 5 Conclusion

This study assessed the performance of CDA, and a discrete variant called DDA, at downscaling the 2D Rayleigh-Bénard convection at various spatial and temporal resolutions of observations. The analysis was conducted in an uncertain setting, where noisy observations of a reference solution are downscaled and statistically analyzed. The CDA and DDA downscaled fields were contrasted qualitatively and quantitatively. In particular, the performances of these algorithms were quantified by systematically analyzing the impact of observation noise level and the temporal and spatial resolutions of the

observational grid, namely by quantifying the skill scores of the downscaled temperature and velocity fields and characterizing their statistics.

The skill scores of the downscaled fields from both CDA and DDA algorithms decay to a plateau, which reflects the convergence of the downscaled solution towards the reference field. Specifically, the CDA skill scores generally decay rapidly and monotonically, as opposed to the DDA skill scores, which exhibit a non-monotonic behavior before reaching a plateau. The plateaus of the DDA skill scores were at lower values compared to the CDA skill scores, meaning that errors of the DDA downscaled fields become smaller than those of CDA fields after the algorithm converges.

A thorough analysis of the functional dependence of the expected $l_2$-error of the downscaled fields ($\Lambda$) was performed on the noise level, temporal and spatial resolutions of the observation grid. $\Lambda$ was observed to vary quadratically with the variance of the observational noise for both CDA and DDA, in agreement with the theoretical results of [22]. Moreover, the CDA results suggest that $\Lambda$ is proportional to the observation's temporal resolution, as opposed to a nonlinear dependence in the case of DDA. Specifically, for small and large temporal resolutions, $\Lambda$ is smaller when downscaling with CDA, however, for intermediate temporal resolutions, DDA downscaled fields result in a smaller $\Lambda$. Finally, our results illustrate that $\Lambda$ varies quadratically with the temporal resolution of the observational data in the case of CDA as opposed to the more sensitive DDA error where $\Lambda$ was observed to vary cubically with the temporal resolution of the observations. The distribution of the downscaled fields were also examined by applying a statistical hypothesis test. Specifically, after reaching the plateau, the downscaled fields pass the Kolmogorov-Smirnov test at almost all tested locations in the domain, suggesting that they are normally distributed.

The precision of the estimated skill scores was also assessed following standard bootstrapping. The results suggest that bootstrapping as few as 30 realizations yields a distribution of the average RRMSE with a standard deviation comparable to that obtained from a 200 realizations. Moreover, the mean downscaled solution achieves a lower error (roughly three times) than any of the downscaled solutions for both CDA and DDA. When the temperature noise level is large, temperature observations degrade the quality of the downscaled fields and discarding them yields a smaller RRMSE. On the other hand, assimilating more accurate temperature observations helps reach the plateau faster and results in more truthful downscaled fields with smaller RRMSE.

Future work will extend the methodology to three-dimensional flows, address the combined effects of model
and observation errors, and assess the performance of both CDA and DDA in comparison to other ensemble-based data assimilation techniques, such as the Ensemble Kalman Filter. We further plan to perform a fundamental analysis of the behavior of the DDA algorithm, namely, to try to establish rigorous error estimates in settings involving observation and model noise.

**Acknowledgments** Research reported in this publication was supported by the Office of Sponsored Research (OSR) at King Abdullah University of Science and Technology (KAUST) CRG Award No. OSR-CRG2020-4336 and the Virtual Red Sea Initiative (Award No. REP/1/3268-01-01).

**Data Availability** The datasets generated during the current study are available from the corresponding author on reasonable request.

## Declarations

**Conflict of Interests** Authors declare that they have no conflicts of interest.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

1. Kim, J., Moin, P., Moser, R.: Turbulence statistics in fully developed channel flow at low reynolds number. J. Fluid Mech. **177**, 133–166 (1987). https://doi.org/10.1017/S0022112087000892

2. Wilby, R.L., Wigley, T.M.L.: Downscaling general circulation model output: a review of methods and limitations. Progress in Physical Geography: Earth and Environment **21**(4), 530–548 (1997). https://doi.org/10.1177/030913339702100403

3. Katavouta, A., Thompson, K.R.: Downscaling ocean conditions with application to the gulf of maine, scotian shelf and adjacent deep ocean. Ocean Model. **104**, 54–72 (2016). https://doi.org/10.1016/j.ocemod.2016.05.007

4. Shin, S.-I., Alexander, M.A.: Dynamical downscaling of future hydrographic changes over the northwest atlantic ocean. J. Clim. **33**(7), 2871–2890 (2020). https://doi.org/10.1175/JCLI-D-19-0483.1

5. Dosio, A., Panitz, H.-J., Schubert-Frisius, M., Lüthi, D.: Dynamical downscaling of cmip5 global circulation models over cortex-africa with cosmo-clm: evaluation over the present climate and analysis of the added value. Climate Dynam. **44**(9), 2637–2661 (2015). https://doi.org/10.1007/s00382-014-2262-x

6. Desamsetti, S., Dasari, H.P., Langodan, S., Titi, E.S., Knio, O., Hoteit, I.: Efficient dynamical downscaling of general circulation models using continuous data assimilation. Q. J.

Roy. Meteorol. Soc. **145**(724), 3175–3194 (2019). https://doi.org/10.1002/qj.3612

7. Huth, R.: Statistical downscaling of daily temperature in central europe. Journal of Climate **15**(13), 1731–1742 (2002). https://doi.org/10.1175/1520-0442(2002)015<1731:SDODTI>2.0.CO;2

8. Tisseuil, C., Vrac, M., Lek, S., Wade, A.J.: Statistical downscaling of river flows. J. Hydrol. **385**(1), 279–291 (2010). https://doi.org/10.1016/j.jhydrol.2010.02.030

9. Jha, S.K., Mariethoz, G., Evans, J., McCabe, M.F., Sharma, A.: A space and time scale-dependent nonlinear geostatistical approach for downscaling daily precipitation and temperature. Water Resour. Res. **51**(8), 6244–6261 (2015). https://doi.org/10.1002/2014WR016729

10. Laflamme, E.M., Linder, E., Pan, Y.: Statistical downscaling of regional climate model output to achieve projections of precipitation extremes. Weather Clim. Extrem. **12**, 15–23 (2016). https://doi.org/10.1016/j.wace.2015.12.001

11. von Storch, H., Langenberg, H., Feser, F.: A spectral nudging technique for dynamical downscaling purposes. Monthly Weather Rev. **128**(10), 3664–3673 (2000). https://doi.org/10.1175/1520-0493(2000)128<3664:ASNTFD>2.0.CO;2

12. Kidson, J.W., Thompson, C.S.: A comparison of statistical and model-based downscaling techniques for estimating local climate variations. J. Clim. **11**(4), 735–753 (1998). https://doi.org/10.1175/1520-0442(1998)011<0735:ACOSAM>2.0.CO;2

13. Altaf, M.U., Titi, E.S., Gebrael, T., Knio, O.M., Zhao, L., McCabe, M.F., Hoteit, I.: Downscaling the 2d bénard convection equations using continuous data assimilation. Comput. Geosci. **21**(3), 393–410 (2017). https://doi.org/10.1007/s10596-017-9619-2

14. Roux, R.L., Katurji, M., Zawar-Reza, P., Quénol, H., Surtan, A.: Comparison of statistical and dynamical downscaling results from the WRF model. Environ Model Softw **100**, 67–73 (2018). https://doi.org/10.1016/j.envsoft.2017.11.002

15. Azouani, A., Olson, E., Titi, E.S.: Continuous data assimilation using general interpolant observables. J. Nonlinear Sci. **24**(2), 277–304 (2013). https://doi.org/10.1007/s00332-013-9189-y

16. Farhat, A., Jolly, M.S., Titi, E.S.: Continuous data assimilation for the 2d bé,nard convection through velocity measurements alone. Physica D: Nonlinear Phenomena **303**, 59–66 (2015). https://doi.org/10.1016/j.physd.2015.03.011

17. Farhat, A., Johnston, H., Jolly, M., Titi, E.S.: Assimilation of nearly turbulent rayleigh-bénard flow through vorticity or local circulation measurements: a computational study. J. Sci. Comput. **77**(3), 1519–1533 (2018). https://doi.org/10.1007/s10915-018-0686-x

18. Farhat, A., Lunasin, E., Titi, E.S.: Data assimilation algorithm for 3d bénard convection in porous media employing only temperature measurements. J. Math. Anal. Appl. **438**(1), 492–506 (2016). https://doi.org/10.1016/j.jmaa.2016.01.072

19. Jolly, M.S., Martinez, V.R., Olson, E.J., Titi, E.S.: Continuous data assimilation with blurred-in-time measurements of the surface quasi-geostrophic equation. Chinese Annals of Mathematics, Series B **40**(5), 721–764 (2019). https://doi.org/10.1007/s11401-019-0158-0

20. Gesho, M., Olson, E., Titi, E.S.: A computational study of a data assimilation algorithm for the two-dimensional navier-stokes equations. Commun. Comput. Phys. **19**(4), 1094–1110 (2016). https://doi.org/10.4208/cicp.060515.161115a

21. Pei, Y.: Continuous data assimilation for the 3d primitive equations of the ocean. Commun. Pure Appl. Anal. **18**(2), 643–661 (2019). https://doi.org/10.3934/cpaa.2019032

22. Bessaih, H., Olson, E., Titi, E.S.: Continuous data assimilation with stochastically noisy data. Nonlinearity **28**(3), 729–753 (2015). https://doi.org/10.1088/0951-7715/28/3/729
23. Hayden, K., Olson, E., Titi, E.S.: Discrete data assimilation in the lorenz and 2d navier–stokes equations. Physica D: Nonlinear Phenomena **240**(18), 1416–1425 (2011). https://doi.org/10.1016/j.physd.2011.04.021

24. Celik, E., Olson, E., Titi, E.S.: Spectral filtering of interpolant observables for a discrete-in-time downscaling data assimilation algorithm. SIAM J. Appl. Dyn. Syst. **18**(2), 1118–1142 (2019). https://doi.org/10.1137/18M1218480

25. Curry, J.H., Herring, J.R., Loncaric, J., Orszag, S.A.: Order and disorder in two- and three-dimensional bénard convection. J. Fluid Mech. **147**, 1–38 (1984). https://doi.org/10.1017/S0022112084001968

26. Paul, M.R., Einarsson, M.I., Fischer, P.F., Cross, M.C.: Extensive chaos in rayleigh-bénard convection. Phys. Rev. E **75**, 045203 (2007). https://doi.org/10.1103/PhysRevE.75.045203

27. Le Quéré, P.: Accurate solutions to the square thermally driven cavity at high rayleigh number. Comput. Fluids **20**(1), 29–41 (1991). https://doi.org/10.1016/0045-7930(91)90025-D

28. Le Maître, O.P., Reagan, M.T., Najm, H.N., Ghanem, R.G., Knio, O.M.: A stochastic projection method for fluid flow: Ii. random process. J. Comput. Phys. **181**(1), 9–44 (2002). https://doi.org/10.1006/jcph.2002.7104

29. Weidauer, T., Pauluis, O., Schumacher, J.: Cloud patterns and mixing properties in shallow moist rayleigh–bénard convection. New J. Phys. **12**(10), 105002 (2010). https://doi.org/10.1088/1367-2630/12/10/105002

30. Chillà, F., Schumacher, J.: New perspectives in turbulent rayleigh-bénard convection. The European Physical Journal E 35(7). https://doi.org/10.1140/epje/i2012-12058-1 (2012)

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.