Data Assimilation in Large-Prandtl Rayleigh-B\'enard Convection from Thermal Measurements

TL;DR

This paper investigates a data assimilation method for large-Prandtl Rayleigh-Bénard convection using temperature observations, providing theoretical guarantees and numerical validation, with implications for geophysical and high-pressure gas applications.

Contribution

It introduces a rigorous framework for data assimilation in large-Prandtl convection, identifying conditions for synchronization and validating them through simulations.

Findings

Synchronization often occurs beyond theoretical conditions

Conditions are practically relevant only at extremely large Prandtl numbers

Numerical results confirm theoretical predictions

Abstract

This work applies a continuous data assimilation scheme---a particular framework for reconciling sparse and potentially noisy observations to a mathematical model---to Rayleigh-B\'enard convection at infinite or large Prandtl numbers using only the temperature field as observables. These Prandtl numbers are applicable to the earth's mantle and to gases under high pressure. We rigorously identify conditions that guarantee synchronization between the observed system and the model, then confirm the applicability of these results via numerical simulations. Our numerical experiments show that the analytically derived conditions for synchronization are far from sharp; that is, synchronization often occurs even when the conditions of our theorems are not met. We also develop estimates on the convergence of an infinite Prandtl model to a large (but finite) Prandtl number generated set of observations. Numerical simulations in this hybrid setting indicate that the mathematically rigorous results are accurate, but of practical interest only for extremely large Prandtl numbers.

Figures

Tables

Table 1. Table 1 : Minimal values of μ 𝜇 \mu that result in synchronization within about 0.005 units of simulation time for the given Ra Ra \operatorname{Ra} with N = 32 𝑁 32 N=32 . These values represent the edge of what works when N = 32 𝑁 32 N=32 : a lower μ 𝜇 \mu may not stimulate convergence, but a larger μ 𝜇 \mu will. See Figure 5 for a quadratic least-squares fit of this data.

$\mathrm{Ra}$	$\mu$
$1.1937766\times {10}^7$	$6,300$
$1.6037187\times {10}^7$	$6,500$
$2.1544346\times {10}^7$	$7,900$
$2.8942661\times {10}^7$	$9,400$
$3.8881551\times {10}^7$	$10,700$
$5.2233450\times {10}^7$	$12,700$
$7.0170382\times {10}^7$	$15,300$
$9.4266845\times {10}^7$	$19,400$
$1.26638017\times {10}^8$	$24,900$
$1.70125427\times {10}^8$	$29,900$

Table 2. Table 2 : Minimal values of N 𝑁 N that result in synchronization within about 0.01 units of simulation time for the given Ra Ra \operatorname{Ra} and μ 𝜇 \mu . In general, as long as an appropriate μ 𝜇 \mu is chosen, N = 16 𝑁 16 N=16 is sufficient for physically relevant values of Ra Ra \operatorname{Ra} .

$\mathrm{Ra}$	$\mu$	$N$
$1.1937766\times {10}^7$	$6,300$	$6$
$1.6037187\times {10}^7$	$6,500$	$8$
$2.1544346\times {10}^7$	$7,900$	$8$
$2.8942661\times {10}^7$	$9,400$	$8$
$3.8881551\times {10}^7$	$10,700$	$10$
$5.2233450\times {10}^7$	$12,700$	$10$
$7.0170382\times {10}^7$	$15,300$	$10$
$9.4266845\times {10}^7$	$19,400$	$10$
$1.26638017\times {10}^8$	$24,900$	$12$
$1.70125427\times {10}^8$	$29,900$	$14$

Table 3. Table 3 : Minimal values of Pr Pr \operatorname{Pr} that result in synchronization within about 0.01 units of simulation time for the given Ra Ra \operatorname{Ra} and μ 𝜇 \mu with N = 32 𝑁 32 N=32 . The relationship here is much less precise than those described in Tables 1 and 2 (for example, Ra ≈ 2.89 × 10 7 Ra 2.89 superscript 10 7 \operatorname{Ra}\approx 2.89\times 10^{7} is a bit of an outlier), but it is interesting to note that lower Ra Ra \operatorname{Ra} seem to require larger Pr Pr \operatorname{Pr} . Whether or not this is an effect of Ra Ra \operatorname{Ra} increasing or μ 𝜇 \mu increasing, however, is unclear.

$\mathrm{Ra}$	$\mu$	$\mathrm{Pr}$
$1.1937766\times {10}^7$	$10,000$	75
$1.6037187\times {10}^7$	$11,000$	56
$2.1544346\times {10}^7$	$12,000$	56
$2.8942661\times {10}^7$	$13,000$	100
$3.8881551\times {10}^7$	$14,000$	42
$5.2233450\times {10}^7$	$18,000$	42
$7.0170382\times {10}^7$	$20,000$	42
$9.4266845\times {10}^7$	$25,000$	31
$1.26638017\times {10}^8$	$32,000$	31
$1.70125427\times {10}^8$	$40,000$	31

Equations

$$\frac{du}{dt}=F(u),$$

$$\frac{dv}{dt}=F(v)-\mu (I_h(v)-I_h(u)),$$

$$v(0)=v_0,$$

$${\partial }_t\mathbf{u}+(\mathbf{u}\cdot \nabla )\mathbf{u}=\nu \Delta \mathbf{u}-\nabla p+\alpha g{\mathbf{e}}_{d}T,\nabla \cdot \mathbf{u}=0,{\partial }_{t}T+(\mathbf{u}\cdot \nabla )T-\kappa \Delta T=0,$$

$$\mathbf{u}{|}_{x_d=0}=\mathbf{u}{|}_{x_d=h}=0,T{|}_{x_d=0}=\delta T,T{|}_{x_d=h}=0,\mathbf{u},T\text{ are L-periodic in }{x}_1,x_{d-1},$$

$$\mathrm{Pr}:=\frac{\nu }{\kappa },\mathrm{Ra}:=\frac{\alpha g(\delta T)h^3}{\nu \kappa }.$$

$$\begin{array}{rl}\frac{1}{\mathrm{Pr}}\left[{\partial }_{t^{\prime }}{\mathbf{u}}^{\prime }+({\mathbf{u}}^{\prime }\cdot {\nabla }^{\prime }){\mathbf{u}}^{\prime }\right]-{\Delta }^{\prime }{\mathbf{u}}^{\prime }& =-{\nabla }^{\prime }{p}^{\prime }+\mathrm{Ra}{\mathbf{e}}_d{T}^{\prime },{\nabla }^{\prime }\cdot {\mathbf{u}}^{\prime }=0,{\mathbf{u}}^{\prime }({\mathbf{x}}^{\prime },0)={\mathbf{u}}_0^{\prime }({\mathbf{x}}^{\prime })\\ {\partial }_{t^{\prime }}{T}^{\prime }+{\mathbf{u}}^{\prime }\cdot {\nabla }^{\prime }{T}^{\prime }-{\Delta }^{\prime }{T}^{\prime }& =0,T^{\prime }({\mathbf{x}}^{\prime },0)=T_0^{\prime }({\mathbf{x}}^{\prime })\end{array}$$

$${\mathbf{u}}^{\prime }{|}_{x_d^{\prime }=0}={\mathbf{u}}^{\prime }{|}_{x_d^{\prime }=1}=0,T^{\prime }{|}_{x_d^{\prime }=0}=1,T^{\prime }{|}_{x_d^{\prime }=1}=0,{\mathbf{u}}^{\prime },T^{\prime }\text{ are L-periodic in }{x}_1^{\prime },x_{d-1}^{\prime },$$

$$\begin{array}{c}-\Delta \mathbf{u}=-\nabla p+\mathrm{Ra}{\mathbf{e}}_{d}T,\nabla \cdot \mathbf{u}=0,\\ {\partial }_{t}T+\mathbf{u}\cdot \nabla T-\Delta T=0,T(\mathbf{x},0)=T_0(\mathbf{x}),\\ \mathbf{u}{|}_{x_d=0}=\mathbf{u}{|}_{x_d=1}=0,T{|}_{x_d=0}=1,T{|}_{x_d=1}=0,\mathbf{u},T\text{ are L-periodic in }{x}_1,x_{d-1}.\end{array}$$

$$\begin{array}{c}\frac{1}{\mathrm{Pr}}[{\partial }_t\tilde{\mathbf{u}}+(\tilde{\mathbf{u}}\cdot \nabla )\tilde{\mathbf{u}}]-\Delta \tilde{\mathbf{u}}=-\nabla \tilde{p}+\mathrm{Ra}{\mathbf{e}}_d\tilde{T},\nabla \cdot \tilde{\mathbf{u}}=0,\tilde{\mathbf{u}}(\mathbf{x},0)={\tilde{\mathbf{u}}}_0(\mathbf{x})\\ {\partial }_t\tilde{T}+\tilde{\mathbf{u}}\cdot \nabla \tilde{T}-\Delta \tilde{T}=-\mu P_N(\tilde{T}-T),\tilde{T}(\mathbf{x},0)={\tilde{T}}_0(\mathbf{x}),\\ \tilde{\mathbf{u}}{|}_{x_d=0}=\tilde{\mathbf{u}}{|}_{x_d=1}=0,\tilde{T}{|}_{x_d=0}=1,\tilde{T}{|}_{x_d=1}=0,\tilde{\mathbf{u}},\tilde{T}\text{ are L-periodic in }{x}_1,x_{d-1},\end{array}$$

$$\begin{array}{c}-\Delta \tilde{\mathbf{u}}=-\nabla \tilde{p}+\mathrm{Ra}{\mathbf{e}}_d\tilde{T},\nabla \cdot \tilde{\mathbf{u}}=0,\\ {\partial }_t\tilde{T}+\tilde{\mathbf{u}}\cdot \nabla \tilde{T}-\Delta \tilde{T}=-\mu P_N(\tilde{T}-T),\tilde{T}(\mathbf{x},0)={\tilde{T}}_0(\mathbf{x}),\\ \tilde{\mathbf{u}}{|}_{x_d=0}=\tilde{\mathbf{u}}{|}_{x_d=1}=0,\tilde{T}{|}_{x_d=0}=1,\tilde{T}{|}_{x_d=1}=0,\tilde{\mathbf{u}},\tilde{T}\text{are L-periodic in }{x}_1,x_{d-1},\end{array}$$

$$\mathcal{F}$$

$${\mathcal{F}}_{\sigma }^d$$

$$H$$

$$V$$

$$W$$

$$⟨\mathbf{u},\mathbf{v}⟩=\sum _{i=1}^3{\int }_{\Omega }{u}_i(\mathbf{x})v_i(\mathbf{x})d\mathbf{x},⟨f,g⟩={\int }_{\Omega }f(\mathbf{x})g(\mathbf{x})d\mathbf{x},$$

$$⟨f,g{⟩}_{H^k}:=\sum _{|\gamma |=0}^k⟨D^{\gamma }f,D^{\gamma }g⟩,$$

$$(f,g{)}_H:=⟨f,g⟩,(f,g{)}_V:=⟨\nabla f,\nabla g⟩,(f,g{)}_W:=(f,b{)}_V+\sum _{|\gamma |=2}⟨D^{\gamma }f,D^{\gamma }g⟩.$$

$$\begin{array}{rl}& W\hookrightarrow V\hookrightarrow H\hookrightarrow H^{\prime }\hookrightarrow V^{\prime }\hookrightarrow W^{\prime },\\ & \mathcal{W}\hookrightarrow \mathcal{V}\hookrightarrow \mathcal{H}\hookrightarrow {\mathcal{H}}^{\prime }\hookrightarrow {\mathcal{V}}^{\prime }\hookrightarrow {\mathcal{W}}^{\prime }.\end{array}$$

$$\{({\lambda }_n,{\varphi }_n(\mathbf{x})){\}}_{n=1}^{\infty }$$

$$f(\mathbf{x},t)=\sum _{n=1}^{\infty }{f}_n(t){\varphi }_n(\mathbf{x}),f_n(t)=⟨f(t),{\varphi }_n⟩={\int }_{\Omega }f(\mathbf{x},t){\varphi }_n(\mathbf{x})d\mathbf{x},$$

$$⟨{\varphi }_i,{\varphi }_j⟩={\delta }_{ij}=\{\begin{array}{ll}1& \text{if }i=j,\\ 0& \text{if }i\ne j.\end{array}$$

$$P_N(f)=\sum _{n=1}^N{f}_n{\varphi }_n,Q_N(f)=[I-P_N](f)=\sum _{n=N+1}^{\infty }{f}_n{\varphi }_n,$$

$$⟨P_N(f),Q_N(f)⟩$$

$$\Vert P_N(f){\Vert }^2+\Vert Q_N(f){\Vert }^2$$

$$-\Delta \mathbf{u}+\nabla p=\mathbf{f},$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\mathbf{u}{|}_{x_3=0}=\mathbf{u}{|}_{x_3=1}=0,\mathbf{u}\text{ is L-periodic in }{x}_1\text{ and }{x}_2,$$

$$\Vert \mathbf{u}{\Vert }_{H^2}+\Vert p{\Vert }_{H^1}\le C\Vert \mathbf{f}\Vert ,$$

$${\psi }^{+}:=\max \{\psi ,0\},{\psi }^{-}:=\max \{-\psi ,0\}.$$

$${\partial }_{t}T+\mathbf{u}\cdot \nabla T-\Delta T=0,T(\mathbf{x},0)=T_0(\mathbf{x}),$$

$$T{|}_{x_d=0}=1,T{|}_{x_d=1}=0,T\text{is a.e. L-periodic in}{x}_1,x_{d-1},$$

$$0\le \bar{T}(t)\le 1,\eta =(T-1{)}^{+}-T^{-},\Vert \eta (t)\Vert \le C_0{e}^{-t}.$$

$$\mathbf{u}\in L^{\infty }(0,\tau ;\mathcal{W}),T\in L^{\infty }(0,\tau ;L^2(\Omega ))\cap L^2(0,\tau ;H^1(\Omega ))\cap C_w([0,\tau ];L^2(\Omega )),$$

Full text

**Keywords: Data Assimilation, Rayleigh-Bénard Convection, Large Prandtl Limit

MSC2010: 76E06, 62M20, 35Q35 **

1. 1 Introduction
2. 2 The Rayleigh-Bénard System and Nudging Equation
3. 3 Mathematical Background
4. 4 Infinite-Prandtl Assimilation
5. 5 Finite-Prandtl Assimilation
6. 6 A Scenario of Model Error
7. 7 Conclusions and Outlook
8. Acknowledgments

1 Introduction

In order to make accurate predictions, numerical models for geophysical processes require establishing accurate initial conditions. Data assimilation is used to estimate weather or ocean (or any other geophysical) variables by incorporating the real world data into the mathematical system to obtain an accurate initialization. One of the classical methods of data assimilation, see, e.g., [HA76, Dal91, LBM91, LSN82, VH89, SS90, ZNLD92, SB93, VLDP03, AB08, BLSZ13, AOT14], is to insert observational measurements directly into a model as the latter is being integrated in time (also known as nudging or newtonian relaxation). There is a significant amount of recent literature concerning the mathematically rigorous analysis of nudging algorithms for data assimilation developed for hydrodynamic equations with a particular focus on weather and climate systems. Recently, a nudging scheme, known as 3DVAR, was studied in [BLSZ13] in the case where observables are given as noisy Fourier modes, and in [AOT14], which successfully accommodates a larger class of observables that, in particular, includes the more physically relevant cases of nodal values and volume elements; see also [BOT15], where observational error is accounted for. In these articles, rigorous proofs are obtained for the synchronization of the approximating signal with the true signal that corresponds to the observations, using the two-dimensional (2D) incompressible Navier-Stokes equations (NSE) as a paradigm.

The data assimilation algorithm analyzed in [BLSZ13, AOT14] can be described as follows: suppose that $u(t)$ represents a solution of some dynamical system governed by an evolution equation of the type

[TABLE]

where the initial state of the system, $u(0)=u_{0}$, is unknown. We would like to accurately track this solution $u(t)$ as $t$ increases notwithstanding our uncertainty in $u_{0}$. Let $I_{h}(u(t))$ represent an interpolant operator based on the observations of the system at a coarse spatial resolution of size $h$, for $t\in[0,T]$. We then construct a solution $v(t)$ from the observations that satisfies the equations

[TABLE]

where $\mu>0$ is a relaxation (nudging) parameter and $v_{0}$ can be prescribed as an arbitrary initial condition. We then take $v(t)$ as prediction of $u(t)$ which we anticipate becomes more accurate as $t$ (and therefore the amount of observed data $I_{h}(u(t))$) increases.

The algorithm designated by (1.2) was designed to work for dissipative dynamical systems of the form (1.1) that are known to have global-in-time solutions, a finite-dimensional global attractor, as well as a finite set of determining parameters (see, e.g., [FP67, FT84, FT91, JT92b, JT92a, CJT97, HT97, FMRT01] and references therein). Typically in these settings, following the ideas in [FP67], lower bounds on $\mu>0$ and upper bounds on $h>0$ can be derived such that the approximate solution $v(t)$ converges to the reference solution $u(t)$ as $t\to\infty$. This was initially demonstrated for the 2D Navier-Stokes equations in [BLSZ13, AOT14].

Numerous further studies, both analytical and numerical, have been carried out for the algorithm (1.2), illustrating its broad scope of applicability. For instance, the nudging approach has been validated for models including the 2D magnetohydrodynamic system [BHLP18], the 2D surface quasi-geostrophic equation [JMT17], three-dimensional (3D) Brinkman-Forchheimer-Extended Darcy model [MTT16], and 3D simplified Bardina model [AB18]. The practically and physically relevant scenarios of discrete-time and time-averaged observables was studied in [FMT16, BMO18, JMOT18]; more recently, it was shown in [BFMT18] that this nudging algorithm is capable of synchronizing the statistics propagated by the flow as they are observed only on a coarse-mesh scale; the efficacy of this algorithm for assimilating actual data sampled from a regional domain encompassing most of Northern Africa and the Middle East was recently tested in [DDL*+*18]. We refer the reader to [Dal91] for a summary on the use of data assimilation in practical forecasting and [Kal03] for a comprehensive text on numerical weather prediction where nudging has been employed.

Regarding related numerical studies, [GOT16] demonstrated in the case of the 2D NSE that the number of observables required for synchronization using (1.2) is much lower in practice than what has been deemed sufficient by the rigorous analysis. In the setting of the 2D RB system, numerical studies were carried out in [ATG*+*17], and then in [FJJT18] for nearly turbulent flows using vorticity and local circulation measurements. We emphasize that the numerical experiments carried out in the present article are for moderately turbulent flows whose dynamics are significantly more complex than the regime of two-cell convection rolls that [ATG*+*17] was restricted to. Moreover our studies are carried out to a similar high degree of numerical precision as found in [FJJT18]. We refer the reader to [ANLT16, FET17, DLMB18, CHL18] for various other studies in the context of turbulent flows such as how one can leverage the nudging scheme to infer unknown parameters of the flow. In [BM17, IMT18, MT18] analytical studies on the various modes of synchronization of the algorithm (1.2) and on certain variants on its numerical discretization were carried out.

The earth system is heated from within and cooled by the atmosphere or ocean at the earth’s surface. On geological time scales the mantle’s motion can be modeled as a fluid. The big difference between the temperature of the top mantle and the bottom mantle is a major source of the convective motion (fluid motion driven by temperature difference). The full compressible, temperature-dependent viscous equations of motion that ostensibly describe flow in the mantle [TS14] are currently beyond the reach of a rigorous mathematical analysis, but a first order approximation to this system is adequately described by taking the infinite Prandtl limit [Wan04a, Wan04b, Wan05, Wan07, Wan08a, Wan08b, FGHR15] of the Rayleigh-Bénard (RB) system first described in [Ray16] by Lord Rayleigh. We recall that the Prandtl number represents the ratio of the kinematic viscosity to the thermal diffusivity. Since the original formulation of the minimal mathematical model in [Ray16], extensive research has sought to quantify its dynamical evolution, [Lor63, Ahl74, AB78, BdBAC91], and resultant large spatial and long temporal scale impact of convective flow, see [AGL09, LX10] for example. Despite the seeming simplicity of the RB system, there remains open questions regarding the exact nature of the convective heat transport, and the impact and nature of boundary layers at physically relevant values (see e.g. [AGL09]). To further complicate matters, mantle convection is far more nuanced than Rayleigh-Bénard convection, having several other unanswered questions, in addition to the well-known open problems in the latter setting. Unlike the low Prandtl number setting, experimental investigations of mantle convection are not practical, so numerical simulations provide one of the only avenues to investigate these issues. Of fundamental concern in such simulations is the dependence of the simulation on the initial condition and/or true physical setting, which would ideally be accompanied by physical observations. The collection of observational data from the mantle is an onerous inverse problem that obfuscates much of the desired resolution in time and space [TS14], both in the sense of physical accessibility, and due to the relevant time scales involved. Indeed, the evolution of the mantle is on millenial time scales. Despite remarkable advances in imaging technologies, observations of the mantle are sparse and prone to noise and are insufficient to determine the mantle’s state. Thus, the development of the advanced real-time prediction systems that are capable of depicting and predicting the mantle’s state is necessary to gain insight into the dynamics in the earth’s interior.

Access to observations from earth’s mantle is limited. Geophysicists have decent observations at the surface (top layer) of plate velocities and heat flow distribution and have probes of mantle temperature where volcanism occurs. Since the motion of the tectonic plates is very slow (the relative movement of the plates typically ranges from zero to 100 mm annually), we can assume that the top plate is stationary and with fixed temperature. Thus, a realistic forecast model for the dynamics of earth’s mantle will employ sparse thermal observations only. In the context of atmospheric and oceanic physics, data assimilation algorithms where some state variable observations are not available as an input, have been studied in [CHJ69, EB87, GSY77, GHA78] for simplified numerical forecast models. Particularly related to the study carried out in this article, it was shown in the case of the 2D NSE [FLT16a] and the 2D RB system in [FLT17] that measurements on just a single component of velocity is sufficient to obtain synchronization. Charney’s question in [CHJ69, GHA78, GSY77] asks whether temperature observations are enough to determine the entire dynamical state of the system. In [GSY77], an analytical argument suggested that Charney’s conjecture is correct, in particular, for a shallow water model. Further numerical testing in [GHA78] affirmed that it is not certain whether assimilation with temperature data alone will yield initial states of arbitrary accuracy. The authors in [ATG*+*17] concluded that assimilation using coarse temperature measurements only will not always recover the true state of the full system. It was observed that the convergence to the true state using temperature measurements only is actually sensitive to the amount of noise in the measured data as well as to the spacing (the sparsity of the collected data) and the time-frequency of such measured temperature data. Rigorous justification for Charney’s Conjecture was provided in [FLT16c] in the case of the 3D Planetary Geostrophic model. Earlier, for the specific setting of 3D convection in a porous medium, where inertial effects can be ignored in the fluid velocity, it was shown in [FLT16b] that temperature measurements alone suffice to determine the velocity field. By comparison, the thrust of the analysis performed here is to establish the conclusions analogous to [FLT16b] while accounting for these inertial effects within the regime of a finite but large Prandtl number.

We consider the nudging approach both analytically and through numerical experiments to explore the range of applicability of the technique in this geophysically interesting context of large Prandtl convective systems. Ultimately, we accomplish the following:

We develop a nudging data assimilation scheme for both large and infinite Prandtl number Rayleigh-Bénard convection in the traditional simplified three-dimensional box geometry (cf. (2.4),(2.5) and (2.9) below) with observations in the temperature field only. In Section 2 we formally introduce the governing equations for the 3D RB system and then in Section  3, we provide the mathematical framework within which our analysis is performed, as well as the relevant well-posedness results for the 3D RB system. We then establish rigorous estimates on the convergence rates for the simpler case of $\Pr=\infty$ in Section 4. The case of large, but finite $\Pr$ is addressed in Section 5. 2. 2.

Perform high-resolution direct numerical simulations (DNS) on the two-dimensional version of this problem for moderately turbulent flows, in an effort to shed some light on the practical applicability of the rigorous estimates. In particular, we probe the values of the relaxation parameter and the number of required modes for the nudging scheme to converge. This is done in Sections 4.2 and 5.2, immediately after their respective mathematical analysis. 3. 3.

Consider a practical scenario of ‘model error’, in which the assimilated variables are nudged “incorrectly.” Specifically, we assume that the modeling system corresponds to the infinite Prandtl system nudged by data corresponding to a finite Prandlt system (2.4), (2.5). This situation is studied both analytically and numerically in Section 6.

We note that the choice of finitely many Fourier modes as the manifestation of our observables is made for ease of both exposition and numerical implementation. We reserve establishing estimates on the convergence rates for more general observables to a subsequent study.

2 The Rayleigh-Bénard System and Nudging Equation

This initial section recalls the Rayleigh-Bénard (RB) system and its non-dimensional formulation. We then present the precise form of the nudging algorithm which we will study in the sequel.

The Rayleigh-Bénard system for convection originates from the Boussinesq equations for an incompressible fluid with appropriate boundary conditions. The Boussinesq system over a $d$-dimensional domain, where $d=2,3$, $\Omega=[0,\tilde{L}]^{d-1}\times[0,h]$, is given by

[TABLE]

where $\mathbf{u}=(u_{1},\dots,u_{d})$ is the velocity vector field, $p$ is the scalar pressure field, and $T$ denotes the temperature of a buoyancy driven fluid. The parameter $\nu$ denotes the kinematic viscosity of the fluid, $\kappa$ its thermal diffusivity, $\alpha$ the thermal expansion coefficient, $g$ denotes the constant gravitational force and $\mathbf{e}_{d}$ is a constant vector anti-parallel to the gravitational force. To model convection, (2.1) is then supplemented by the boundary conditions

[TABLE]

where $\delta T$ is a fixed constant determined by the (relative) strength of the bottom heating. The relations (2.1), (2.2) together constitute the 3D RB system for convection. We note that variations on these boundary conditions are applicable to the earth’s mantle, but as our focus is on the convergence of the data assimilation scheme, we assume that such variations have secondary effects.

Non-dimensionalized variables

As is customary, we work with non-dimensionalized variables. The system (2.1) is rescaled using $h$ as a length scale, $\delta T$ as the temperature scale, and the diffusive scale $\frac{\kappa}{h^{2}}$ as the time scale. The relevant non-dimensional physical parameters for the system are the Prandtl number, $\operatorname{Pr}$, and Rayleigh number, $\operatorname{Ra}$, which are defined as

[TABLE]

This leads to non-dimensionalized variables over the rescaled domain $\Omega^{\prime}=[0,L]^{d-1}\times[0,1]$, $d=2,3$, which satisfy

[TABLE]

with the boundary conditions

[TABLE]

For notational simplicity, we will drop the ′ in all that follows.

As previously alluded to, the physical setting of interest in this article is the earth’s mantle where the Prandtl number is large, namely, on the order of $10^{25}$. Upon formally setting $\operatorname{Pr}=\infty$ in the system (2.4)–(2.5), one arrives at

[TABLE]

We note that the initial velocity $\mathbf{u}(\mathbf{x},0)$ is then determined by $T_{0}$ and the corresponding momentum equations. Although there are several additional physical effects relevant to mantle convection that are omitted from the Rayleigh-Bénard model considered here, we consider (2.9) to be an appropriate “zeroth order” representation of mantle convection; it provides the starting point and test model for mantle convection simulations (cf. [BBC*+*89]). Although we anticipate eventually extending the results of the current investigation to more realistic models of mantle convection, we believe a more in-depth understanding of the problem at hand is a necessary first step.

Nudging setup

The idea, following [BLSZ13, AOT14] and sketched in the introduction, is to nudge the assimilated system as in (1.2) with a projection of the ‘truth’ that represents the exactly realizable observations of the original system. More precisely, the nudging is accomplished by introducing an affine feedback control term to the original ‘forecast’ model (1.1), whose purpose is to enforce the asymptotic convergence of the solution of the assimilated system (1.2) towards that of the original system (1.1), but only on the scales at which the observations are made; it is this ‘relaxed’ imposition that ensures the practicality of the nudging scheme. For this article, our ‘truth’ is assumed to be represented by (2.4), (2.5) or by (2.9).

Let $(\mathbf{u},T)$ satisfy (2.4), (2.5) over $\Omega=[0,L]^{2}\times[0,1]$, from which we have obtained partial observations in the form of finitely many Fourier coefficients corresponding to wave-numbers $|k|\leq N$, for some integer $N>0$. Let $(\tilde{\mathbf{u}},\tilde{T})$ denote the assimilated or modeled system variables, which satisfies

[TABLE]

where $P_{N}$ denotes the projection onto Fourier wave-numbers $|k|\leq N$ (see (3.8) below). Its infinite Prandtl counterpart is given by

[TABLE]

where $(\mathbf{u},T)$ comes from either (2.4), (2.5) or (2.9).

One of the basic goals of this paper is to show that $\tilde{T}-T,\tilde{p}-p\rightarrow 0$ and $\tilde{\mathbf{u}}-\mathbf{u}\rightarrow\mathbf{0}$ as $t\rightarrow\infty$ in the appropriate space for specific conditions on $\mu$ and the number of projected modes $N$ relative to $\operatorname{Ra}$ and $\operatorname{Pr}$. This indicates that for specified Rayleigh $\operatorname{Ra}$ and Prandtl $\operatorname{Pr}$ numbers one can determine a sufficiently large number of modes and a sufficiently large nudging parameter $\mu$ to ensure that the assimilated system $(\tilde{\mathbf{u}},\tilde{T})$ will asymptotically match the true system.

3 Mathematical Background

For the sake of completeness, this section presents some preliminary material and notation commonly used in the mathematical study of hydrodynamic systems, in particular in the study of the NSE and the Euler equations for incompressible fluids. For more detailed discussion on these topics, we refer the reader to, e.g., [CF88, Tem97].

Let $\Omega=[0,L]^{d-1}\times[0,1]$, where $d=2,3$ and we denote the spatial variable $\mathbf{x}=(x_{1},\ldots,x_{d})$. We consider the function spaces

[TABLE]

where $H^{1}(\Omega)$ and $H^{2}(\Omega)$ denote the classical Sobolev class of first-order and second-order weakly differentiable functions over $\Omega$, respectively. We use the notation $X^{\times d}$ to denote the $d$-fold product of a set $X$, and $\langle\cdot,\cdot\rangle$ to denote the usual $L^{2}$ inner product over $\Omega$,

[TABLE]

for ${\bf u}=(u_{1},\dots,u_{d}),{\bf v}=(v_{1},\dots,v_{d})\in{\mathcal{H}}$, and $f,g\in L^{2}(\Omega)$. The inner product on $H^{k}(\Omega)$, $k=1,2,\dots$, is given by

[TABLE]

where $\gamma=(\gamma_{1},\dots,\gamma_{d})$ is a multi-index and $D^{\gamma}=(\partial_{x_{1}}^{\gamma_{1}},\dots,\partial_{x_{d}}^{\gamma_{d}})$. The spaces $H,V,W$ are then endowed with a Hilbert space structure, whose respective inner products are given by

[TABLE]

The spaces $\mathcal{H},\mathcal{V},\mathcal{W}$ have analogous Hilbert space structures. We denote by $H^{\prime},V^{\prime},W^{\prime}$ and $\mathcal{H}^{\prime},\mathcal{V}^{\prime},\mathcal{W}^{\prime}$ the dual spaces of $H,V,W$, and $\mathcal{H},\mathcal{V},\mathcal{W}$, respectively. We then have the following continuous injections:

[TABLE]

In what follows we will denote the $L^{2}(\Omega)$ norm by $\|\cdot\|$. For all other Banach spaces $X$, e.g. $L^{p}(\Omega)$, for $p\not=2$, $H^{k}(\Omega)$ etc., we denote the associated norms explicitly as $\|\cdot\|_{X}$.

Let

[TABLE]

denote the orthonormal eigenpairs corresponding to the Laplace operator $-\Delta$ on the domain $\Omega$ supplemented with the mixed horizontally periodic-vertically Dirichlet boundary condition as in (3.1). Then each $f\in H^{2}(\Omega)$ can be expressed in terms of the eigenfunctions as

[TABLE]

and the eigenfunctions satisfy the orthogonality relation

[TABLE]

For each $N\geq 0$, define the projections $P_{N},Q_{N}$ by

[TABLE]

where $I$ is the identity operator. In other words, $P_{N}$ is a truncation of the eigenfunction expansion, and $Q_{N}$ is its orthogonal complement. The orthogonality of the eigenbasis yields the identities

[TABLE]

for any $f\in H^{2}(\Omega)$.

We next recall the following well-known a priori estimates for Stokes’ equations

[TABLE]

cf. [CF88] and the drift-diffusion equation, (3.15) where the advecting velocity field is divergence free as in e.g. [Tem97, Wan05]. We adapt these results here to our notation and setting as follows

Lemma 3.1**.**

Let $d=2,3$ and $\mathbf{f}\in H^{\times d}$. There exists a unique $\mathbf{u}\in\mathcal{W}$ and (up to constants) $p\in V$ such that (3.11) is satisfied. Moreover, there exists a constant $C=C(\Omega)>0$ such that

[TABLE]

The next lemma follows from the theory of linear transport equations in [DL89] and a variant of the maximum principle proved in [FMT87] .We will refer to the following notation for the “positive” and “negative parts” of a function:

[TABLE]

Lemma 3.2**.**

Let $d=2,3$ and $\tau>0$. Let $\mathbf{u}\in L^{1}(0,\tau;\mathcal{V})$ and $T_{0}\in L^{\infty}(\Omega)$ be a.e. $L$-periodic in $x_{1},x_{d-1}$ and $T_{0}|_{x_{d}=0}=0,T_{0}|_{x_{d}=1}=1$ (in the sense of trace). Suppose that $T\in L^{\infty}(0,\tau;L^{\infty}(\Omega))\cap L^{2}(0,\tau;H^{1}(\Omega))$ satisfies

[TABLE]

where the boundary values on $\{x_{d}=0\}\cup\{x_{d}=1\}$ are interpreted in the sense of trace. Then there exists a constant $C_{0}=C_{0}(\Omega,\lVert T_{0}\rVert)>0$ and functions $\bar{T},\eta$ such that $T=\bar{T}+\eta$ and satisfy

[TABLE]

for all $t>0$.

For the system (2.9), when $\operatorname{Pr}=\infty$, the velocity field, $\mathbf{u}$, is determined by the evolution of $T$. The well-posedness of this system then follows in a standard way, for either dimension $d=2,3$, and its solution satisfies the estimates stated in Lemma 3.1 and 3.2. We formally state this result as the following theorem.

Theorem 3.1**.**

Let $d=2,3$. Suppose that $T_{0}\in L^{\infty}(\Omega)$ is a.e. in $\Omega$, that $T_{0}$ is a.e. $L$-periodic in $x_{1},x_{d-1}$ and $T_{0}|_{x_{d}=0}=0,T_{0}|_{x_{d}=1}=1$ (in the sense of trace). Then there exists a unique $(\mathbf{u},T)$ satisfying (2.9) such that

[TABLE]

for all $\tau>0$. Moreover, there exist positive constants $\gamma_{0}=\gamma_{0}(\Omega,\lVert T_{0}\rVert)$ and $C_{0}=C_{0}(\Omega,\|T_{0}\|)$, such that

[TABLE]

for all $t\geq 0$, where $\gamma_{0}^{\prime}:=|\Omega|^{1/2}+C_{0}e^{-t}$.

We consider a change of variable, denoted by $\theta({\bf x},t)$, that represents the fluctuation of the temperature around the steady state background temperature profile $1-x_{d}$:

[TABLE]

The functional setting determined by (3.2)-(3.4) accommodates a rigorous mathematical analysis for the perturbed variable $\theta$. The results derived for $\theta$ are then transferred naturally to the desired results for the original variable $T$. We appeal to [Wan05, Wan07] for the global existence and eventual regularity of suitable weak solutions for (2.4)-(2.5), as well as the existence of the “global attractor” for the dynamics, although we will not make explicit use of this fact in this article. We will also say that a solution $(\mathbf{u},T)$ of (2.4), (2.5) is regular on $[0,\tau]$ if $(\mathbf{u},\theta)\in L^{\infty}(0,\tau;\mathcal{V})\times L^{\infty}(0,\tau;V)$. If $\tau=\infty$, we say that the solution is a global regular solution.

Theorem 3.2** ([Wan05, Wan07]).**

Let $d=2,3$. Recalling the notation (3.2)–(3.4) let $(\mathbf{u}_{0},\theta_{0})\in\mathcal{H}\times H$ and $T_{0}:=\theta_{0}+(1-x_{d})$.

1. (i)

(Global Existence of Weak Solutions) For any $\tau>0$, there exists $(\mathbf{u},T)$ such

[TABLE]

where $\theta$ is related to $T$ by the relation (3.16). This $(\mathbf{u},T)$ satisfies (2.4)–(2.5) in the usual weak sense, and maintains the initial condition $(\mathbf{u}_{0},T_{0})$. Moreover, the following energy inequality and maximum principle are satisfied for all $t\leq\tau$:

[TABLE]

where $(T-1)^{+},T^{-}$ are defined as in (3.14). 2. (ii)

(Eventual Regularity of Weak Solutions) There exists a constant $K_{0}>0$ such that if $\operatorname{Pr}\operatorname{Ra}^{-1}\geq K_{0}$, then there exists a time $\tau^{*}>0$ for which all suitable weak solutions corresponding to $(\mathbf{u}_{0},\theta_{0})\in\mathcal{H}\times H$ become regular solutions on $[\tau^{*},\infty)$. In particular, when $d=3$, there exist constants $\kappa_{1},\kappa_{2},\kappa_{3}>0$ and $\kappa_{0}^{\prime\prime},\kappa_{1}^{\prime\prime},\kappa_{2}^{\prime\prime}>0$, depending on $\Omega$, such that

[TABLE]

for all $t\geq\tau^{*}$.

For the remainder of the article, we will therefore assume that $(\mathbf{u},T)$, $(\tilde{\mathbf{u}},\tilde{T})$ have evolved for a sufficiently long time, so that $(\mathbf{u},T)$ is a regular solution to either (2.4)–(2.5) or (2.9). Physically speaking, the set-up of our study assumes that reality is represented exactly by a solution to (2.4)–(2.5) or (2.9) and that we have been observing the system after the point in time at which it has become globally regular. Thus, for the purposes of our analysis, we will henceforth make the following standing hypotheses for the remainder of the paper.

Standing Hypotheses.

Let $d=2,3$ and $\operatorname{Ra}\geq 1$. Let $\gamma_{0}>0$ be the constants from Theorem 3.1. When $\Pr=\infty$, we assume:

1. $(I1)$

$T_{0}$ a.e. periodic in $x_{1},x_{d-1}$ and $T_{0}|_{x_{d}=0}=0$, $T_{0}|_{x_{d}=1}=1$ (in the sense of trace); 2. $(I2)$

$T_{0}\in L^{\infty}(\Omega)$; 3. $(I3)$

$(\mathbf{u},T)$ is the unique global solution to (2.9) corresponding to $T_{0}$ and guaranteed by Theorem 3.1. In particular, $(\mathbf{u},T)$ satisfies

[TABLE]

for all $t>0$.∎

On the other hand, let $K_{0}>0$, $\kappa_{1},\kappa_{2},\kappa_{3}>0$, and $\kappa_{0}^{\prime\prime},\kappa_{1}^{\prime\prime},\kappa_{2}^{\prime\prime}>0$ be the constants in Theorem 3.2 $(ii)$. When $\Pr<\infty$, we assume:

1. $(F1)$

$\Pr\operatorname{Ra}^{-1}\geq K_{0}$; 2. $(F2)$

$(\mathbf{u},T)$ is the unique regular solution to (2.4), (2.5); 3. $(F3)$

When $d=3$, $\mathbf{u}(t)$ satisfies

[TABLE]

for all $t\geq 0$. 4. $(F4)$

When $d=3$, $T$ satisfies

[TABLE]

for all $t\geq 0$, where $\kappa_{j}^{\prime}=\kappa_{j}^{\prime\prime}+2|\Omega|^{1/2}+1$, $j=0,1,2$.∎

Remark 3.1**.**

Since we have non-dimensionalized our variables, we point out that although the bounds in $(F3),(F4)$ are derived for the case $d=3$, they are also valid for the case $d=2$, up to constants, provided that $\operatorname{Ra}\geq 1$, which have have assumed as one of our standing hypotheses.

4 Infinite-Prandtl Assimilation

We will first treat the Data Assimilation problem for the infinite Prandtl number system (2.9), (2.17). Hence, throughout this section we assume that $\operatorname{Pr}=\infty$, i.e., $(I1)-(I3)$ holds.

A rigorous mathematical analysis is performed in dimensions $d=2,3$, while the numerical component of our studies are carried out for $d=2$. Due to the structure of the nudged system (2.17), we do not have a maximum principle. Instead, we require only that (2.17) have a well-defined solution in the weak sense, which one can do by establishing that the differences $\mathbf{w}=\tilde{\mathbf{u}}-\mathbf{u}$, $S=\tilde{T}-T$ satisfy their respective evolution in the weak sense. Since one of the relevant a priori estimates to this end are performed below for the proof of Theorem 4.2, we simply state this result as the following theorem.

Theorem 4.1**.**

Let $\mu>0$ and $\{\lambda_{n}\}_{n=1}^{\infty}$ be as in (3.7). Let $\tilde{T}_{0}\in L^{2}(\Omega)$ a.e. in $\Omega$ such that ${\tilde{T}}_{0}$ is a.e. $L$-periodic in $x_{1},x_{d-1}$ with ${\tilde{T}}_{0}|_{x_{d}=0}=0,{\tilde{T}}_{0}|_{x_{d}=1}=1$ (in the sense of trace). Suppose $N>0$ satisfies $\frac{1}{4}\lambda_{N}\geq\mu$. Then there exists a unique $(\tilde{\mathbf{u}},{\tilde{T}})$ satisfying (2.17) in the weak sense such that

[TABLE]

*for all $\tau>0$. *

4.1 Synchronization

We are now ready to state and prove the main theorem of this section. We show that (2.9) and (2.17) synchronize under certain conditions detailed below.

Theorem 4.2** (Infinite-Prandtl Synchronization).**

Let $\mu>0$ and $\{\lambda_{n}\}_{n=1}^{\infty}$ be as in (3.7). Let $N>0$ satisfy $\frac{1}{4}\lambda_{N}\geq\mu$ and ${\tilde{T}}_{0}$ be given as in Theorem 4.1. Let $(\tilde{\mathbf{u}},\tilde{T})$ be the corresponding unique solution to (2.17) guaranteed by Theorem 4.1. There exists a constant $C_{0}=C_{0}(\Omega,\|T_{0}\|)>0$ such that if

[TABLE]

then for all $t\geq 0$

[TABLE]

for $C_{1}=\|\tilde{T}_{0}-{T}_{0}\|$.

Proof.

Let $S=\tilde{T}-T$, $\mathbf{w}=\tilde{\mathbf{u}}-\mathbf{u}$, and $q=\tilde{p}-p$. Subtracting (2.9) from (2.17) yields the system

[TABLE]

with the initial condition $S(\mathbf{x},0)=\tilde{T}_{0}(\mathbf{x})-T_{0}(\mathbf{x})$. The momentum equation in (4.6) satisfies Lemma 3.1, so there exists a constant $C>0$ such that

[TABLE]

Therefore, to establish (4.2), it is sufficient to show that $S\rightarrow 0$ with an exponential rate in $L^{2}(\Omega)$.

Upon multiplying the $S$ equation in (4.6) by $S$ and integrating over $\Omega$, we obtain

[TABLE]

Assume that $\mu$ is chosen sufficiently large so that (4.1) holds, where $C_{0}>0$ is, as of yet, unspecified. After integrating by parts, an application of the Sobolev embedding $H^{2}(\Omega)\hookrightarrow L^{\infty}(\Omega)$, (4.7), and then $(I3)$ of the Standing Hypotheses imply

[TABLE]

On the other hand, due to (3.10) and the inverse Poincaré inequality, and since $N>0$ satisfies $\frac{1}{4}\lambda_{N}\geq\mu$, it follows that

[TABLE]

Joining (4.8) with (4.2) and (4.2) and combining like terms, we arrive at

[TABLE]

Finally, by Gronwall’s inequality and the second condition in (4.1), we deduce that

[TABLE]

which completes the proof upon choosing $C_{0}=C\gamma_{0}^{\prime}$. ∎

Theorem 4.2 shows that given $\operatorname{Ra}$ and reasonable boundary conditions $T_{0}$ and $\tilde{T}_{0}$, we can always choose $\mu$ and $N$ large enough so that $(\mathbf{u},T)$ and $(\tilde{\mathbf{u}},\tilde{T})$ eventually match in the infinite-time limit. In the next subsection, we use numerical simulations to verify that this is indeed the case.

4.2 Numerical Results at Infinite Prandtl number

Rather than focusing on a handful of highly turbulent three-dimensional simulations, for the same computational cost we consider more detailed and extensive simulations in the two-dimensional setting where $\Omega=[0,L]\times[0,1]$ with coordinates $\mathbf{x}=(x_{1},x_{3})$ in order to search through the relevant parameters. We simulate (2.9) and (2.17) using a stream function formulation with Dedalus [BVO*+*16], a Python package that uses pseudospectral methods to solve partial differential equations on spectrally representable domains. All of the simulations in this section and all that follow are completed with a 4-stage 3rd order Runge-Kutta implicit-explicit time stepping scheme that treats the linear terms implicitly and the nonlinear terms explicitly. Each simulation runs with $L=4$ and $256$ Fourier grid points in the horizontal and $128$ Chebyshev points in the vertical with a standard $3/2$ dealiasing. All simulations were run until the time-averaged Nusselt number and other pertinent statistics were temporally well-converged for all cases considered here.

Choosing the initial conditions $T_{0}$ and $\tilde{T}_{0}$ in our numerical experiments is nontrivial. For $T_{0}$, we have two options.

• •

Set $T_{0}(x_{1},x_{3})=1-x_{3}+\varepsilon(x_{1},x_{3})$ for some small perturbation function $\varepsilon:\Psi\rightarrow\mathbb{R}$. This begins the simulation close to the conductive state $1-x_{3}$, which—though not physically relevant—is a suitable starting point for initial experiments.

• •

Load $T_{0}$ from the final state of a previous simulation with similar parameters. If this previous simulation has run long enough, $T$ will thus begin in a reasonable state for the new set of parameters.

Finding a suitable choice for $\tilde{T}_{0}$ is even less obvious. Setting $\tilde{T}_{0}(x_{1},x_{3})=1-x_{3}$ assumes no prior knowledge of $T$, and results in an overly stiff initialization due to the strength of the initial temperature difference, and is therefore numerically impractical. In addition, initiating $T$ at this conductive state ignores observations taken at the initial time $t=0$ which would not be advantageous in practice. On the other hand, setting $\tilde{T}_{0}=P_{N}(T_{0})$ for large $N$ results in a very weak temperature difference so that (2.9) and (2.17) are nearly synchronized at $t=0$. As a compromise, we set $\tilde{T}_{0}$ as a low-mode projection of $T_{0}$ such as $\tilde{T}_{0}=P_{4}(T_{0})$. This permits some initial knowledge of the true state without driving it directly to the ‘truth’ immediately.

Values of the Rayleigh number $\operatorname{Ra}$ that are of interest for mantle convection are typically between $10^{7}$ and $10^{8}$ [TS14, BBC*+*89]. However, the greater value for $\operatorname{Ra}$, the more computationally expensive the simulation, and it is numerically unstable to start $T$ in a state such as $T_{0}(x_{1},x_{3})\approx 1-x_{3}$ at large $\operatorname{Ra}$. We therefore select logarithmically spaced $\operatorname{Ra}$ values from $10^{4}$ to $10^{9}$ and run a simulation for each $\operatorname{Ra}$, one after another. For the very first simulation ($\operatorname{Ra}=10^{4}$), we set $T_{0}(x_{1},x_{3})=1-x_{3}+\varepsilon(x_{1},x_{3})$, as discussed above, and evolve the systems forward until reaching a nearly steady state. Thereafter, we set $T_{0}$ as the final $T$ from the previous simulation and increment increase $\operatorname{Ra}$. This yields realistic initial conditions for any given $\operatorname{Ra}$ in the range considered. Figure 1 shows an example of $T_{0}$ and $\tilde{T}_{0}$, together with subsequent end states of $T$ and $\tilde{T}$ for $\operatorname{Ra}\approx 5.22\times 10^{7}$.

To measure the synchronization of (2.9) and (2.17), we keep track of the $L^{2}(\Omega)$, $H^{1}(\Omega)$, and $H^{2}(\Omega)$ norms of $\tilde{T}-T$ and $\tilde{\mathbf{u}}-\mathbf{u}$ throughout the simulation. Each error is normalized by dividing by the norms of the truth system’s variables. For example, to measure the temperature difference in $L^{2}(\Omega)$, we compute $\frac{\|(\tilde{T}-T)(t)\|}{\|T(t)\|}$ at each simulation time $t$. For simplicity, we write these errors without the denominators in all that follows.

Figure 2 shows how the tracked error norms decrease in the simulation from Figure 1. Though the analysis indicates that these norms should converge to $\varepsilon_{\text{machine}}\approx 10^{-16}$, seeing the errors decrease to about $10^{-9}$ is numerically satisfactory in this setting, given the format in which the norms are calculated, i.e. some errors do accumulate in the calculation of the norm itself.

4.2.1 Dependencies of nudging parameter

Previous studies, [ATG*+*17, FJJT18], of data assimilation for Rayleigh-Bénard convection fixed the value of the nudging parameter at $\mu=1$ in their experiments and instead focused on the effects of other parameters, such as the number of projected modes. However, Theorem 4.2 requires that $\mu$ be proportional to $\operatorname{Ra}^{2}$, so we do not expect synchronization with $\mu=1$ as $\operatorname{Ra}$ increases. Indeed, Figure 3 shows that for $\operatorname{Ra}$ in the range of $4\times 10^{7}$ to $2\times 10^{8}$, using $\mu=1$ results in either extremely slow convergence, or in no convergence at all.

To explore the relationship between $\operatorname{Ra}$, $\mu$, and $N$ needed for synchronization to occur, we fix two of the parameters at a time and run several simulations with various values for the remaining parameter. First, we fix $\operatorname{Ra}$, set $N=32$, and vary $\mu$ until synchronization can be observed within $0.005$ units of simulation time (a simulation of this length usually requires several thousand iterations). Table 1 records the smallest value of $\mu$ where synchronization was observed; Figure 4 shows the convergence rates for a few different $\mu$ values when $\operatorname{Ra}\approx 3.89\times 10^{7}$ and $\operatorname{Ra}\approx 7.02\times 10^{7}$.

Given the requirement (4.1) from Theorem 4.2, it is surprising to discover that—based on Table 1—the relationship between $\mu$ and $\operatorname{Ra}$ is more linear than quadratic. Indeed, a least squares fit of the data to a general quadratic of the form $f(x)=a+bx+cx^{2}$ yields coefficients of $a\approx 4.03\times 10^{3}$, $b\approx 1.77\times 10^{-4}$, and $c\approx-1.37\times 10^{-13}$ (essentially $c=0$). See Figure 5.

4.2.2 Relation between relaxation and number of observables

The relaxation parameter $\mu$ is a system parameter without a clear physical interpretation. The number of projected Fourier modes $N$, on the other hand, indicates the amount of data that is “visible” to the assimilating system. Therefore, a lower bound on $N$ represents how much data is required in order to maintain an accurate model. Fixing $\mu$ as given by Table 1, we decrease $N$ to see how it affects synchronization. See Figure 6. Note that as expected, the less modes retained in the projection, the slower synchronization is achieved, and if a sufficiently low number of modes ($N\leq 6$) is observed then synchronization never occurs. In Table 2 we summarize our findings concerning the number of modes $N$ needed as a function of $\operatorname{Ra}$ (given a suitably large choice for $\mu$).

4.2.3 Summary

The numerical simulations verify that over physically relevant values of $\operatorname{Ra}$, we can pick $\mu$ large enough to nudge (2.9) toward (2.17); in particular, the experiments confirm Theorem 4.2 and show that the conditions stated therein are stronger than necessary. Furthermore, synchronization can be obtained with a relatively low number of modes even for large $\operatorname{Ra}$, as long as $\mu$ is chosen large enough. We took the approach of varying $\mu$ first and then $N$, but it is just as feasible to vary $\mu$ for a fixed $\operatorname{Ra}$ and $N$.

5 Finite-Prandtl Assimilation

We now turn to the case when $\operatorname{Pr}$ is finite, but large. Thus, throughout this section, we will assume that $\operatorname{Pr}<\infty$, so that we automatically have that $(F1)-(F4)$ holds.

As in Section 4, we first state the relevant well-posedness result for the nudged equation. We then rigorously establish synchronization of the nudged signal with the true signal, then proceed with a numerical study for the two-dimensional setting.

It is not immediately clear how to properly adapt the argument for (2.9)–(2.17) in [Tem97, Wan05] to establish a maximum principle (see Lemma 3.2) for the assimilating variable $\tilde{T}$ due to the presence of the additional term $-\mu P_{N}(\tilde{T}-T)$ in (2.17) (cf. [JMST18] for a case where it is necessary to develop a maximum principle in spite of this term). However, since we are assuming $(\mathbf{u},T)$ is a regular solution to (2.4)–(2.5), it is straightforward to verify global existence and uniqueness for the associated nudged system by instead considering the corresponding system for the difference $(\mathbf{w},S)$, where $\mathbf{w}=\tilde{\mathbf{u}}-\mathbf{u}$ and $S=\tilde{T}-T$. In this setting, we need only appeal to a maximum principle for $T$, rather than for $\tilde{T}$. The well-posedness of the nudged system then follows in a standard fashion, so that we opt to only state the result here and refer the reader to [FJT15] for the appropriate details.

Theorem 5.1**.**

Let $\mu>0$ and $\{\lambda_{n}\}_{n=1}^{\infty}$ be as in (3.7). Let $(\tilde{\mathbf{u}}_{0},\tilde{T}_{0}-(1-x_{d}))\in\mathcal{V}\times V$. Suppose $N>0$ satisfies $\frac{1}{4}\lambda_{N}\geq\mu$. Then there exists a unique solution $(\tilde{\mathbf{u}},\tilde{T})$ satisfying (2.13) such that

[TABLE]

for all $\tau\geq 0$.

5.1 Synchronization

We will establish the finite Prandtl analog to Theorem 4.2. For this, it will be convenient to establish the following stability estimate first.

Lemma 5.1**.**

Assuming the conditions of Theorem 5.1, there exists a constant $C_{0}=C_{0}(\Omega)>0$ such that if

[TABLE]

then

[TABLE]

holds for $t\geq 0$.

Proof.

Similar to the analysis performed for Theorem 4.2, we consider the solution $(\mathbf{w},S)$ to the difference system

[TABLE]

Then, multiplying the second equation by $S$ integrating and arguing as in (4.2) we obtain

[TABLE]

On the other hand, an analogous energy calculuation for $\mathbf{w}$ yields

[TABLE]

We estimate the right hand side of (5.7) using Hölder’s inequality, the Gagliardo-Nirenberg interpolation inequality, and Young’s inequality to obtain

[TABLE]

Next, for $II$ in (5.8), we have from Hölder’s inequality, the Sobolev embedding $L^{6}(\Omega)\hookrightarrow H^{1}$, Gagliardo-Nirenber interpolation inequality, and Young’s inequality that

[TABLE]

For $III$, we use the Cauchy-Schwarz inequality to obtain

[TABLE]

Combining the bounds for $I$-$III$, we have

[TABLE]

Thus, by the second condition in (5.1), it follows that

[TABLE]

Therefore, by Gronwall’s inequality the desired bound (5.2) now follows. ∎

Theorem 5.2**.**

Assume the conditions of Lemma 5.1. There exists a constant $C_{1}=C_{1}(\Omega)>0$ such that if

[TABLE]

then

[TABLE]

*holds for all $t\geq 0$, where $C_{2}=\|\tilde{\mathbf{u}}_{0}-\mathbf{u}_{0}\|+\|\tilde{T}_{0}-T_{0}\|$. *

Proof.

By $(F3)$ of the Standing Hypotheses, it follows that

[TABLE]

By $(F4)$ and the Sobolev embedding, $H^{1}(\Omega)\hookrightarrow L^{6}(\Omega)$, it follows that

[TABLE]

Now upon combining (5.12), (5.13), it follows that there exists $C=C(\Omega)>0$ such that

[TABLE]

Then, having chosen $\operatorname{Pr}$ according to (5.10), upon taking square roots, we deduce (5.11) from and Lemma 5.1, as desired. ∎

Remark 5.1**.**

Observe that (5.10) requires the Prandtl number to be very large to ensure that synchronization occurs. Indeed, for $\operatorname{Ra}\sim 10^{7}$ this condition indicates that $\Pr\gtrsim 10^{35}$ which would place the Prandtl number in a regime beyond the one that occurs for the earth’s mantle. In addition, the lower bound on $\mu$ stated here would imply that $\mu\gtrsim\Pr\operatorname{Ra}^{2}$ which would yield a very stiff problem numerically, particularly if $\Pr\sim\operatorname{Ra}^{5}$ as indicated. Gratefully, as seen in the next subsection, these rigorous estimates are pessimistic, and synchronization is achieved for much lower values of $Pr$ and $\mu$ than are indicated here.

5.2 Numerical Results

The numerical results are similar to those presented in Section 4.2, i.e. the same spatial discretization and time-stepping algorithm are employed here, but now we also consider variations in $\operatorname{Pr}$. In addition, for finite $\operatorname{Pr}$ the velocity field is no longer slaved directly to the temperature field, requiring us to specify an initial condition for the velocity as well as the temperature. This is done as in the infinite $\operatorname{Pr}$ case by setting the initial state of the assimilated system to a low-order projection of the ‘true’ system. Simulations at higher Rayleigh numbers are initiated with flow fields from previous simulations at incrementally lower $\operatorname{Ra}$. To ensure that the algorithm is working, we first check that synchronization still occurs for reasonable choices of $\operatorname{Ra}$, $\mu$, and $N$ given a finite $\operatorname{Pr}$, say $\operatorname{Pr}=100$. See Figure 7.

To simplify our exploration, for each value of $\operatorname{Ra}$ listed in Tables 1 and 2, we pick $\mu$ so that the systems synchronize quickly. Then, with $N=32$, we run simulations for logarithmically spaced values of $\operatorname{Pr}$ from $1$ to $100$. Even with these larger-than-necessary choices of $\mu$, convergence is lost for small enough $\operatorname{Pr}$. See Figure 8 for a few additional examples and Table 3 for the chosen $\mu$ and the lowest $\operatorname{Pr}$ where synchronization still occurs.

There are two main points to take away from these results.

• •

Finite-but-large Prandtl data assimilation through only temperature measurements is possible, though it is slower than the infinite Prandtl setting explored in Section 4.

• •

The relationship between $\operatorname{Ra}$, $\mu$, $N$, and $\operatorname{Pr}$ remains unclear and would require further measurements to precisely quantify. However, the fact that the assimilation works at all with reasonably small $\operatorname{Pr}$ indicates that the inequality constraints in Theorem 5.2 are strongly overstated, or that the constants in (5.10) are extremely small.

• •

While the estimates obtained in Theorem 5.2 are clearly pessimistic it is interesting to note that the qualitative behavior is as expected. That is, the finite Prandtl system will synchronize so long as $\operatorname{Pr}$ is sufficiently large. This is not unexpected as the limit of $\operatorname{Pr}\rightarrow 0$ will be dominated by inertial effects in which the temperature and velocity field are nearly decoupled so we would anticipate that temperature observations alone will not suffice to recover the full flow field.

6 A Scenario of Model Error

Finally, we address a realistic scenario of assimilating observables that correspond to solutions of the finite, but large Prandtl number system (2.4), (2.5) into the “incorrect,” albeit more computationally tractable, infinite Prandtl data assimilation system (2.17). Thus, we suppose $(\mathbf{u},T)$ satisfy $(F1)-(F4)$ and simply choose a suitable $\tilde{T}(\mathbf{x},0)=\tilde{T}_{0}(\mathbf{x})$ for (2.17) since the corresponding initial velocity is enslaved by the temperature evolution for the later nudging equation. Before we perform the analysis for the error estimates, we state the well-posedness result corresponding to the nudged equation, whose proof follows along similar lines to that of Theorem 4.1.

Theorem 6.1**.**

Let $\mu>0$, and $\{\lambda_{n}\}_{n=1}^{\infty}$ be as in (3.7). Let $(\mathbf{u},T)$ of (2.4)-(2.5) satisfying $(F1)-(F4)$. Let $\tilde{T}_{0}\in L^{2}(\Omega)$ a.e. in $\Omega$ such that ${\tilde{T}}_{0}$ is a.e. $L$-periodic in $x_{1},x_{d-1}$ with ${\tilde{T}}_{0}|_{x_{d}=0}=0,{\tilde{T}}_{0}|_{x_{d}=1}=1$ (in the sense of trace). Suppose $N>0$ satisfies $\frac{1}{4}\lambda_{N}\geq\mu$. Then there exists a unique $(\tilde{\mathbf{u}},{\tilde{T}})$ satisfying (2.17) in the weak sense such that

[TABLE]

*for all $\tau>0$. *

6.1 Error estimates

Since our data assimilation equation (2.17) does not correspond to the true evolution of the observables (2.4), (2.5), we do not expect to obtain an exact synchronization. Instead, we derive estimates that quantify the maximal error possible, and which will vanish as the Prandtl number is taken increasingly large.

Theorem 6.2**.**

Let $N,\mu>0$ satisfy $\frac{1}{4}\lambda_{N}\geq\mu$ and ${\tilde{T}}_{0}$ be given as in Theorem 6.1. Let $(\tilde{\mathbf{u}},\tilde{T})$ be the corresponding unique solution to (2.17) guaranteed by Theorem 6.1. There exists a constant $C_{0}=C_{0}(\Omega)>0$ such that if

[TABLE]

then there exist positive constants $C_{1}=\|\tilde{T}_{0}-T_{0}\|,C_{2}=C_{2}(\Omega),C_{3}=C_{3}(\Omega)$ such that

[TABLE]

for all $t\geq 0$.

Proof.

Let $S=\tilde{T}-T$, $\mathbf{w}=\tilde{\mathbf{u}}-\mathbf{u}$, $q={\tilde{p}}-p$. Then

[TABLE]

Upon taking the $L^{2}$-inner product of $\mathbf{w},S$ with (6.3), (6.4), respectively, then adding the consequent relations, we obtain

[TABLE]

Upon applying Hölder’s inequality, the Poincaré inequality, the Sobolev embedding theorem, Ladyzhenskaya’s inequality, and Theorem 3.2 $(ii)$, we derive

[TABLE]

Upon combining estimates for $I-IV$ with the conditions in (6.1) and the fact that $\frac{1}{4}\lambda_{N}\geq\mu$, it follows that

[TABLE]

Hence, by Gronwall’s inequality, we arrive at

[TABLE]

Lastly, by Lemma 3.1 and $(F3)$ of the Standing Hypotheses, we have

[TABLE]

We take the square root of (6.5) and add the result to (6.6) to complete the proof.

∎

6.2 Numerical Results

To numerically verify Theorem 6.2 in a way that is consistent with the numerical simulations corresponding to Theorems 4.2 and 5.2, we compare $(\tilde{\mathbf{u}},\tilde{T})$ to $(\mathbf{u},T)$. To begin, consider a simulation with $\operatorname{Ra}\approx 5.22\times 10^{7}$, $\mu=18,000$ $N=32$, and $\operatorname{Pr}=100$. For the finite Prandtl model in Section 5, this set of parameters results in synchronization (see Figure 7). In this situation, however, the synchronization appears to be limited by the $O(\Pr^{-1})$ error from the $\Pr=\infty$ model to the $\Pr<\infty$ reality.

This apparent lack of convergence is expected, however, since Theorem 6.2 only guarantees that the error between $(\tilde{\mathbf{u}},\tilde{T})$ and $(\mathbf{u},T)$ decreases to $O(\operatorname{Pr}^{-1})$ as time increases. Using the same set of parameters, but with larger and larger $\operatorname{Pr}$, results in tighter and tighter synchronization. To more carefully match the results to the statement of Theorem 6.2, we calculate $\operatorname{Ra}^{-1}\|(\tilde{\mathbf{u}}-\mathbf{u})(t)\|_{H^{2}}$ + $\|(\tilde{T}-T)(t)\|$ at each simulation time $t$. See Figure 10.

While Figure 10 is encouraging, it also highlights a weakness in the data assimilation scheme unique to this hybrid setting. In both the previous settings considered here, the rigorous estimates were pessimistic particularly for the large but finite Prandtl case. It appears that the estimates provided here for this hybrid setting are far more true to practice, i.e. while the error is not linear in $\Pr$ as seen in Figure 10 it is certainly dominated by the effects of the Prandtl number. This suggests that the practical success of these types of data assimilation schemes is highly dependent on the data coming from the same model as the simulated system, an unrealistically stringent restriction. Recall that the Boussinesq approximation is effectively a “zeroth order” approximation for the mantle, meaning the true physical system has several complicated secondary effects (some of which are unknown) that are not included in the model. The lack of numerical synchronization in such a simple setting suggests that data assimilation may not be as adaptable to settings where the exact model is not known.

7 Conclusions and Outlook

In Section 4, we examined a data assimilation scheme for the Rayleigh-Bénard system with $\operatorname{Pr}=\infty$ and showed rigorously that synchronization occurs between the data and assimilating equations under certain conditions on the relaxation parameter $\mu$ and the number of projected modes $N$ relative to the Rayleigh number $\operatorname{Ra}$ when measurements of the temperature only are observed. That is, as long as there is enough data (i.e., $N$ is not too small), $\mu$ can be chosen large enough to guarantee synchronization. Though this is a satisfying theoretical result, the numerical experiments in Section 4.2 show that synchronization often occurs under much weaker conditions on $\mu$ and $N$ than Theorem 4.2 requires. In particular, the inequality conditions on $\mu$ is shown to be at least an order of magnitude away from being sharp. In addition, the numerical results also demonstrate situations in which synchronization fails, namely when $\mu$ and/or $N$ are not large enough.

Section 5 shows that synchronization in the temperature measurements only and at finite $\operatorname{Pr}$ is also possible, although the rate of convergence is slower than with $\operatorname{Pr}=\infty$ and the relationship between $\operatorname{Ra}$, $\operatorname{Pr}$, $\mu$, and $N$ needed to achieve synchronization remains somewhat ambiguous. As in the $\operatorname{Pr}=\infty$ case, the conditions imposed on $\mu$ appear quite pessimistic when compared to numerical experiments.

Finally, when the true values are taken from $\operatorname{Pr}<\infty$ simulations, but the assimilating equations use $\operatorname{Pr}=\infty$, the synchronization is highly dependent on $\operatorname{Pr}$, as predicted by the rigorous bounds. This hybrid setting illustrates that the difference between the two systems is dominating the error, rather than the dynamical error in the synchronization process. Although the numerical simulations agree well with the rigorous predictions in this setting, they do indicate a pessimistic outlook for additional settings wherein the exact evolution of the dynamics for a data assimilation system of this type is unknown. In particular as noted above, we have omitted several details in our model of mantle convection that play a vital role in the evolution and may have an effect similar to the difference between finite and infinite $\Pr$. To investigate this further, data assimilation applied to the internally heated convective setting (see [Gol16, WD11, WD12] for example), and possibly the anelastic or compressible convective systems [KLvK*+*10] will be explored.

The current consideration of the difference between the infinite and near-infinite Prandtl number convective systems lends itself to further investigations wherein the assimilating model is different from the physical system wherefrom the observations are obtained. For example, one might consider the effects of imprecisely defined boundary conditions, i.e. what if the observations were obtained from a convective simulation in which the velocity satisfied a Navier-slip condition, but the nudged system was modeled with a no-slip condition? Other variations in the model itself might include slight variations in the geometry between the two systems, and additional terms in the equations themselves such as internal heating mentioned above. The rub of the matter is that data assimilation techniques, if they are meant to apply to physical settings such as weather, climate, and investigations of the earth’s mantle, must consider the fallibility of the model they are relying on, that is do variations in the underlying model itself allow for synchronization of the model with the observed truth?

Acknowledgments

This paper was initiated by conversations between AF and NEGH (and eventually JPW) which took place during a workshop at the Institute for Pure and Applied Mathematics (IPAM), for which all the authors are indebted to IPAM, the National Science Foundation (NSF), and organizers of the program on Mathematics of Turbulence in Fall 2014. NEGH would like to acknowledge the grants NSF-DMS-1816551, NSF-DMS-1313272 and Simons Foundation 515990 which supported this work. SM and JPW acknowledge computational resources and support from the Fulton Supercomputing Laboratory at Brigham Young University where all the simulations were performed, as well as extensive feedback from the dedalus developers (see http://dedalus-project.org/). JPW also acknowledges the hospitality and generosity of the Mathematics Department at Tulane University where some of the finishing touches on this work were carried out.