Jean-Christophe Loiseau
email: loiseau.jc@gmail.com

J.-Ch Loiseau

Data-driven modeling of the chaotic thermal convection in an annular thermosyphon

Keywords: Reduced-order model, System Identification, Chaos, Natural convection

Identifying accurate and yet interpretable low-order models from data has gained a renewed interest over the past decade. In the present work, we illustrate how the combined use of dimensionality reduction and sparse system identification techniques allows us to obtain an accurate model of the chaotic thermal convection in a two-dimensional annular thermosyphon. Taking as guidelines the derivation of the Lorenz system, the chaotic thermal convection dynamics simulated using a high-fidelity computational fluid dynamics solver are first embedded into a low-dimensional space using dynamic mode decomposition. After having reviewed the physical properties the reduced-order model should exhibit, the latter is identified using SINDy, an increasingly popular and flexible framework for the identification of nonlinear continuous-time dynamical systems from data. The identified model closely resembles the canonical Lorenz system, having the same structure and exhibiting the same physical properties. It moreover accurately predicts a bifurcation of the high-dimensional system (corresponding to the onset of steady convection cells) occuring at a much lower Rayleigh number than the one considered in this study.

Introduction

Fluid flows are characterized by high-dimensional nonlinear dynamics that gives rise to rich structures. Despite this apparent complexity, the dynamics often evolves on a lowdimensional attractor defined by a few dominant coherent structures that contain significant energy or are useful for control [START_REF] Holmes | Turbulence, Coherent Structures, Dynamical Systems and Symmetry[END_REF]. Given this property, one might aim to derive or identify reduced-order models that reproduce qualitatively and quantitatively the dynamics of the full system. Over the past decades, identifying robust, accurate and efficient reducedorder models has thus become a central challenge in fluid dynamics and closed-loop flow control [START_REF] Fabbiane | [END_REF]3,4,5,6]. Amidst the numerous flows exhibiting such properties, the buoyancydriven flow in a thermosyphon is of particular interest both from a theoretical and practical point of view. From the engineering perspective, natural convection in closed loops can play an important role in the design of thermal energy systems such as solar heating systems and nuclear reactors. Natural convection is also of crucial importance in numerous geophysical situations such as the mesoscale convective thunderstorms, land and sea breezes resulting from a differential heating between landmass and an adjacent body of water, or Hadley cells in the Earth's atmosphere. From a theoretical point of view finally, natural convection has attracted a lot of attention ever since the seminal work of Edward Lorenz [7] on the derivation of a low-order model for the Rayleigh-Bénard convection.

As in [7], many traditional model reduction techniques are analytical. They rely on prior knowledge of the Navier-Stokes equations and project them onto the span of an othonormal basis of modes, resulting in a dynamical system in terms of the coefficients of this expansion basis. These modes may come from a classical expansion, such as Fourier modes or Tchebyshev polynomials, or they may be data-driven as in the proper orthogonal decomposition (POD) [8,9]. Although such approaches have proven to be quite successful for linear systems [5], they have been applied only with limited success to obtain low-order approximations of nonlinear systems, mostly on flow oscillators [10]. Lately, data-driven approaches have been becoming increasingly popular. They encompass a wide variety of different approaches such the eigenrealization algorithm [11], the dynamic mode decomposition [12,13,14] and its variants, or NARMAX [START_REF] Billings | Nonlinear System Identification: NARMAX Methods in the Time, Frequency, and Spatio-Temporal Domains[END_REF]. Advances in machine learning are also greatly expanding our ability to extract governing dynamics purely from data. One can cite from instance the genetic programming-based system identification proposed by Bongard & Lipson [START_REF] Bongard | [END_REF] and Schmidt & Lipson [17] or the more recent framework based on sparsitypromoting regression techniques, SINDy, proposed by Brunton, Proctor & Kutz [18].

Considering a two-dimensional annular thermosyphon, the aim of the present work is to illustrate how these recent advances in dimensionality reduction and sparsity-promoting regression techniques can be harvested to enable the identification of a low-order Lorenzlike system able to accurately reproduce qualitatively and quantitatively the key features of chaotic natural convection in an annular thermosyphon purely from data. The present paper is organized as follow. First, the exact flow configuration considered is presented in § 2 along with a brief overview of its chaotic dynamics. Then, § 4 briefly introduces the reader to the dimensionality reduction and system identification techniques used, namely dynamic mode decomposition [12,13] and the sparse identification of nonlinear dynamics [18]. Herein, DMD enables us to obtain a low-dimensional representation of the flow's dynamics while SINDy is used to identify an interpretable low-order model correctly approximating these chaotic dynamics. Finally, the key results of this study are presented in § 5 and § 6 while § 7 discusses the new perspectives opened by the present work.

2 Flow configuration, dynamics and statistical analysis Section §2.1 introduces the flow configuration and the non-dimensional numbers characterizing our setup. An overview of the flow dynamics for the set of parameters considered is presented in §2.2. This section also shows strong qualitative evidences that the chaotic dynamics of the flow are closely related to the Lorenz system. Finally, §2.3 summarizes the results of a statistical analysis of the flow dynamics.

Flow configuration and non-dimensional numbers

The geometry considered consists in two concentric circular enclosures, the inner radius being R 1 while the outer one is R 2 as shown in figure 1. In the rest of this work, the ratio of the outer to inner radius is set to

R 2 R 1 = 2.
A constant temperature T 0 is set at the upper walls while the lower ones are set at a temperature T 1 = T 0 + ∆ T , with ∆ T > 0. Hereafter, we will work with a non-dimensionalized temperature θ (x, y) defined as θ (x, y) = T (x, y) -T 0 ∆ T , where x and y are the horizontal and vertical coordinates. Using this non-dimensionalization, the temperature at the lower walls is thus θ w (y < 0) = 1 while the temperature at the upper ones is set to θ w (y > 0) = 0. The origin of our reference frame is chosen as the center of the thermosyphon. Gravity is acting in the vertical direction along -e e e y and is characterized by the gravitational acceleration g.

Assuming the working fluid enclosed within the thermosyphon is Newtonian, it is characterized by its density ρ, its dynamic viscosity µ, its thermal expansion coefficient β and its thermal diffusivity α. Using ∆ T as the temperature scale and R 2 -R 1 as the length scale, we can then define two non-dimensional parameters, namely the Rayleigh number

Ra = ρgβ ∆ T (R 2 -R 1 ) 3 µα
and the Prandtl number

Pr = ν α ,
where ν = µ /ρ is the kinematic viscosity of the working fluid. Throughout this manuscript, the Rayleigh number is fixed to Ra = 17 000 while the Prandtl number is set at Pr = 5.

For the sake of simplicity, we will assume that the flow within the thermosyphon is two-dimensional and incompressible so that the effect of density variations due to temperature can be modeled using the Boussinesq approximation. Under these assumptions, the dynamics of the flow are governed by the following Navier-Stokes equations

∂ u u u ∂t + ∇ • (u u u ⊗ u u u) = -∇p + Pr ∇ 2 u u u + Ra Pr θ e e e y ∂ θ ∂t + (u u u • ∇) θ = ∇ 2 θ ∇ • u u u = 0
where u u u(x x x,t) is the velocity field, p(x x x,t) is the pressure field and θ (x x x,t) is the temperature one. Expressing the equations in this form implies that time has been non-dimensionalized with respect to the diffusive temperature time scale. The mesh consists in 32 spectral elements uniformly distributed in the azimuthal direction and 8 elements uniformly distributed in the radial one. Within each element, Lagrange interpolants of order 7 based on the Gauss-Lobatto-Legendre quadrature points are used in each direction resulting in 16 384 grid points for the velocity and temperature fields. Lagrange interpolants of order 5 based on Gauss-Lobatto quadrature points are used for the pressure. Temporal integration is performed using a third-order accurate scheme and the time-step has been chosen as to satisfy the condition CFL ≤ 0.5 during the whole simulation.

Flow dynamics

The Navier-Stokes equations have been integrated forward in time until a statistical steady state is reached. Figure 2 depicts the instantaneous temperature (first column), radial velocity (middle column) and azimuthal velocity (right columns) fields at three time instants. These instantaneous fields highlight that the flow is made of a single convection cell. This convection cell can either rotate clockwise or counter-clockwise 1 . The dynamics of these seemingly random switches between clockwise and counter-clockwise rotations of the convection cell can be visualized in figure 3(a) depicting a typical time-series of the flow rate m(t) through a cross-section of the thermosyphon. Although the convection cell appears to rotate most of the time in the clockwise direction (m(t) > 0) over the time span of figure 3(a), the empirical probability density function of the flow rate computed over more than 650 diffusive time-units and depicted in figure 3(b) shows that both rotations are equally likely. Figure 3(a) shows moreover that the convection cell appears to oscillate at the same characteristic frequency whether it rotates clockwise or counter-clockwise. Looking at the power spectral density of the time-derivative ṁ(t) of the flow rate in figure 3(c), this characteristic frequency is estimated to be ω 45.36 corresponding to a period τ = 0.138. Although no linear stability analysis of the Navier-Stokes equations has been conducted during this study, we hypothesize that this characteristic frequency is close to the frequency of the leading eigenmode of the linearly unstable steady clockwise (or counter-clockwise) convection cell. Confirmation of this hypothesis is left for future studies as it would require the introduction of additional analyses beyond the scope of the present work.

A natural question to ask is whether these switches between clockwise and counterclockwise rotations of the convection cell are governed by a random stochastic process or by an unknown deterministic one. To answer to this question, one can use tools from dynamical systems theory such as the Takens embedding theorem [19] or the Broomhead-King attractor reconstruction technique [20] which combines ideas from the Takens embedding theorem with the singular value decomposition of a Hankel matrix constructed from timelagged versions of a single time-series. In this work, the latter has been preferred as it enables the reconstruction of the underlying attractor from a single time-series while also providing a simple diagnostic to determine its dimension by looking at the singular values distribution of the Hankel matrix. For more details about this attractor reconstruction technique, interested readers are refered to the original paper [20] as well as to Appendix B. Using the time-series of the flow rate, the attractor reconstructed using the Broomhead-King technique is shown in figure 3 (d). The distribution of singular values (not presented) strongly suggest that this attractor is only three-dimensional. Additionally, the Poincaré section (not shown) highlights that, even though this attractor is embedded within a three-dimensional space, it is almost two-dimensional. Closer inspection of the Poincaré section indicates that it moreover has a fractal structure with a Haussdorf dimension slightly larger than 2, albeit longer simulations would be required to precisely determine this fractal dimension. These empirical results thus strongly suggest that, despite the high-dimensionality of the discretized Navier-Stokes equations, the dynamics of the convection cell are governed by an underlying low-dimensional deterministic process exhibiting chaos. The vast majority of this work is thus dedicated to identifiying this low-dimensional chaotic system from data obtained by direct numerical simulations. It must be noted finally that the attractor reconstructed in figure 3(d) looks strikingly similar to the one obtained when applying the same technique to the famous Lorenz system, see figure 4. This striking resemblance let us think that the unknown low-dimensional system we look for may be very similar to the Lorenz system. This plausible analogy, already noted by various other authors [21,22,23] when studying similar flows, will guide all of the coming analyses.

Statistical analysis

Before diving in the identification of the unknown low-dimensional system governing the dynamics of the convection cell, let us first gain more insights by characterizing more quantitatively the statistics of the system. Figure 5 shows the mean state and the diagonal entries of the corresponding Reynolds stress tensor. These have been obtained by averaging over more than 650 diffusive time-units the various fields once the simulation had reached a statistical steady state. All of these quantities are either symmetric or anti-symmetric with respect to both the vertical and the horizontal axes. Given that the fixed point of the Navier-Stokes equations has not been computed, we cannot conclude regarding the similiarties of the base flow and the mean flow albeit we expect these two to be relatively similar. Looking at the bottom row of figure 5, it is worthy to note moreover that, although the fluctuations of azimuthal velocity are relatively invariant in the azimuthal direction, the temperature fluctuations are essentially localized in the vicinity of y = 0, i.e. where the temperature applied onto the walls jumps from θ w = 1 to θ w = 0.

Let us now look more closely at the organization of the flow when it rotates either in the clockwise or counter-clockwise direction. Figure 6 depicts the the mean flow conditionned on the negative flow rate, i.e.

ū u u CCW (x x x) = E [u u u(x x x,t)|m(t) < 0] ,
As expected, it describes a convection cell rotating in the counter-clockwise direction. Similarly, ū u u CW (x x x) = E [u u u(x x x,t)|m(t) > 0] corresponds to a convection cell rotating in the clockwise direction (not shown). Even though their amplitudes change, the spatial distribution of the Reynolds stresses associated to this conditionally averaged mean flow appears quite similar that of the true mean flow (not shown). These two conditional mean flows are moreover symmetric of one another with respect to the vertical axis, consistent with the observation made from the empirical probability density function of the flow rate that the clockwise and counter-clockwise rotations of the convection cell are equally likely. Within the embedded phase space, these two conditional mean states are located approximately at the center of the empty regions on each side of the attractor (not shown), thus suggesting that they may be very similar to the linearly unstable steady solutions of the Navier-Stokes equations corresponding to the two configurations. From a dynamical system point of view, we expect these two points to be saddle-foci, i.e. linearly unstable fixed points having a single complex conjugate pair of unstable eigenvalues while all the others are stable (i.e. negative real part). Once again, confirmation of this hypothesis is left for future studies.

Finally, let us look back once more at the time-series of the flow rate m(t) and inspect further its statistical properties. For that purpose, we will use symbolic dynamics by intro- Broomhead-King reconstruction of the attractor Fig. 4 (a) Time-series of the x variable of the Lorenz system for (σ , ρ, β ) = (10, 28, 8 /3). The whole timeseries extends from t = 0 to t = 200. Only a subset is shown. (b) Broomhead-King reconstruction of the underlying strange attractor from the time-series shown in (a). The temporal window used for this reconstruction is set to τ = 0.032 non-dimensional time-units corresponding to 32 measurements of the x variable. Each panel of figure (b) shows the projection of the strange attractor onto the planes (v 1 , v 2 ), (v 1 , v 3 ) and (v 2 , v 3 ), respectively, where v i is the i th Broomhead-King coordinate.

x(t) (a) Time-series of x(t) v 1 v 2 v 1 v 3 v 2 v 3 (b)
ducing the telegraph signal s(t). It is defined as

s(t) = sign (m(t)) ,
i.e. s(t) = 1 if the convection cell rotates clockwise and s(t) = -1 otherwise. A time-series of this telegraph signal s(t) along with the corresponding evolution of the flow rate m(t) is shown in figure 7(a). By discarding the oscillatory behaviour in the vicinity of the saddlefoci, the introduction of the telegraph signal s(t) thus enables us to focus our attention onto the switching dynamics only. Looking at the autocorrelation functions of the signals depicted in figure 7(b), it can be observed that the autocorrelation decay observed at short time (i.e. τ < 0.25) is primarily dictated by the switching dynamics. One can additionally look at the statistical distribution of the duration of the clockwise and counter-clockwise rotations. This distribution is shown in figure 7(c). It is close to being a discrete distribution with the peaks being located at constant intervals of time from one another. The most likely duration of rotation in one direction is τ 0 0.12, consistent with the frequency measured from the power spectral density of ṁ(t) in figure 3(c) and the characteristic decay time scale of the autocorrelation functions. This characteristic duration corresponds to the typical time it takes the system to perform a complete oscillation around one of the two sadde-foci. The other peaks in the distribution in figure 7(c) thus correspond to events during which the convection cell oscillates two or more times before it switches direction. The amplitude of the peaks at τ 0 , 2τ 0 , 3τ 0 , etc appears to decay exponentially fast thus highlighting that the convection cell is unlikely to oscillate more than a handful of times before switching direction. All of the statistical properties just described are thus properties that the low-dimensional system we aim to identify would need to verify before we can claim that it captures the essence of the thermosyphon's chaotic dynamics. It should be noted finally that very similar results have been reported in [24] when studying the Lorenz system, hence further highlighting the possible connections that exist between the two systems.

3 From Rayleigh-Bénard to the Lorenz system This section essentially is of pedagogical interest. Given the similarities highlighted in the previous section between our flow configuration and the classical Lorenz system, its aims is to show how one can derive the Lorenz system from the equations governing the Rayleigh-Bénard convection between two infinite horizontal plates. Additionally, this derivation will provides us with guidelines for the identification of the unknown low-dimensional model in §5. Readers already familiar with this derivation may skip directly to §3.3 where the properties of the Lorenz are discussed.

Rayleigh-Bénard convection

The classical setup for Rayleigh-Bénard convection is that of a horizontal layer of fluid heated from below and cooled from above in which the fluid develops a regular pattern of convection cells known as Bénard cells. From a mathematical point of view, the dynamics of the flow are governed by the same set of equations as the thermosyphon with the exception that the length scale used to define the non-dimensional numbers is now chosen as the height H between the two plates. As discussed in the introduction, this classical setup for Rayleigh-Bénard convection serves as a simplified model in numerous fields (e.g. atmospheric convection in meteorology).

Assuming that the lower wall is at temperature θ (y = 0) = 1 and the upper one θ (y = 0) = 0, it can easily be shown that the equations admit the following fixed point

u u u(x x x) = 0 and Θ (x x x) = 1 -y.
This solution corresponds to a pure conducting state.

In order to simplify our derivation of the Lorenz system, let us now introduce the streamfunction ψ(x x x,t) defined such that u u u(x x x,t) = ∇ × ψ(x x x,t)e e e z with appropriate boundary conditions. Rewriting the Navier-Stokes equations in terms of the fluctuations with respect to this base state yields

∂ ∂t ∇ 2 ψ = -J (ψ, ∇ 2 ψ) + RaPr ∂ θ ∂ x + Pr∇ 4 ψ ∂ θ ∂t = -J (ψ, θ ) + ∇ 2 θ + ∂ ψ ∂ x (1)
where the nonlinear terms are expressed as a Jacobian operator given by

J ( f , g) = ∂ f ∂ x ∂ g ∂ y - ∂ g ∂ x ∂ f ∂ y .
Expressing everything in terms of the fluctuations to the base state enables us to simply impose a zero Dirichlet boundary conditions for the temperature fluctuation θ (x x x) on both walls while, as Lorenz, we will impose free-slip boundary conditions for the velocity. Finally, periodic boundary conditions are assumed in the horizontal direction.

Truncated Galerkin expansion and the Lorenz system

Saltzman [25] has shown that ψ(x x x,t) and θ (x x x,t) in Eq. ( 1) can be expressed as infinite Fourier series in both x and y. Rather than using these infinite Fourier series which would yield intractable calculations, let us rather consider a three-terms truncated Galerkin expansion given by ψ(x, y,t) = a(t) sin(πy) sin(kπx)

+ • • • θ (x, y,t) = b(t) sin(πy) cos(kπx) + c(t) sin(2πy) + • • •
The scalars a(t) and b(t) are the amplitudes of the velocity and temperature fields of the convection cells with wavenumber k in the horizontal direction while c(t) describes the modification of the mean temperature profile due to the onset of convection.

Introducing this ansatz into Eq. ( 1) and performing a Galerkin projection yields

da dt = -Pr(k 2 + π 2 )a + PrRa k k 2 + π 2 b db dt = πkac + ka -(k 2 + π 2 )b dc dt = k π 2 ab -4π 2 Pr c.
After some change of variables, the above model can be transformed into the canonical Lorenz system dx dt = σ (yx)

dy dt = x (ρ -z) -y dz dt = xy -β z
where σ plays the role of the Prandtl number, ρ that of the Rayleigh number and β characterizes the aspect ratio of the convection cells. For more details about this derivation, interested readers are referred to [25].

Properties of the Lorenz system

A very large number of studies have been dedicated to the Lorenz system over the past 50 years. The discussion that comes thus does not aim to delve extensively into all of the properties of the Lorenz system, but rather at highlighting the key properties we expect the to-be-identified low-order model to exhibit if the chaotic dynamics of the thermosyphon are indeed Lorenz-like. The first property of interest is that the Lorenz system belong to the class of dissipative dynamical systems, i.e. volumes in phase space contract under the flow f f f (x, y, z). Denoting by V 0 a volume of nearby initial conditions, it can easily be shown that its evolution is governed by

dV dt = -(σ + 1 + β )V
where

-(σ + 1 + β ) = ∇ • f f f (x, y, z).
The solution of this ordinary differential equation is

V (t) = e -(σ +1+β )t V 0 .
Given that σ , β > 0, volumes in phase space thus shrink exponentially fast to a set of measure zero (i.e. fixed points, periodic orbits or strange attractor). When identifying the loworder model in §5, such a condition will impose some constraints on the set of admissible coefficients of the model. The second property of interest is that the Lorenz equations are equivariant with respect to the transformation (x, y, z) → (-x, -y, z).

This equivariance property implies that all solutions are either symmetric or have a symmetric partner. From a physical point of view, it describes the fact that the clockwise and counterclockwise rotations of the convection cell are symmetric from one another and equally likely.

Additionally, as will be discussed in §5, this property also limits the pool of functions one can use for system identification. Finally, the last property of interest is related to the nature of the z-axis. Setting

x = y = 0 ∀t, the Lorenz system reduces to dz dt = -β z and its solution is simply z(t) = e -βt z 0 .

The z-axis thus forms an invariant set and belongs to the stable manifold of the fixed point located at the origin of the phase space. Physically, this mathematical property describes the fact that, starting from a temperature distribution uniform in the horizontal direction and forcing the fluid to be at rest, the system will naturally relaxes exponentially rapidly to the pure conducting state characterized by a linear temperature profile in the vertical direction and uniform in the horizontal one. From our point of view, the three aforementioned mathematical properties are the fundamental ones the to-be-identified system needs to verify (along with the statistical properties presented in §2.3) to be physically consistent. After this pedagogical interlude, the rest of this work now focuses on the numerical methods and the actual identification of our unknown low-order model.

Numerical methods for low-dimensional system identification

In the derivation of his famous equations, Lorenz has used prior knowledge of the governing equations and of their mathematical properties to find a natural system of coordinates (i.e. the leading Fourier-Fourier modes) to represent the state of the system and to derive the equations governing its dynamics within this low-dimensional space. In the present work, the singularity of our temperature boundary condition at y = 0 prevents us for doing so. As a consequence, we thus need to turn our attention to dimensionality reduction and system identification techniques to find a proper low-dimensional embedding of the state and approximate/identify its governing equations within this low-dimensional space. To do so, §4.1 introduces the reader to the dynamic mode decomposition (DMD), an increasingly popular dimensionality reduction proposed by Schmid [13] in 2010. Similarly, §4.2 provides a rapid introduction to the sparse identification of nonlinear dynamics (SINDy), a new framework for the identification of continuous-time systems proposed by Brunton et al. [18]. Using DMD on our data will thus provides us with a low-dimensional representation of the flow's evolution from which a low-order model will be fitted using SINDy.

Dynamic Mode Decomposition

Dynamic mode decomposition (DMD) is an increasingly popular dimensionality reduction technique that has originated from the field of fluid dynamics [12,13]. Since its introduction, numerous variants have been proposed, see for instance [12,13,14,26,27,[START_REF] Kutz | Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems[END_REF][START_REF] Dawson | [END_REF]30,31,32]. It has also been related to Koopman theory [33]. Considering a zero-mean sequence of evenly sampled m-dimensional snapshots {q q q(x x x,t k )} k=1,n , DMD aims at finding a linear operator A A A ∈ R m×m solution to minimize

A A A n-1 ∑ k=1 q q q k+1 -A A Aq q q k 2 2
where q q q k = q q q(x x x,t k ). Introducing the data matrix X X X = q q q 1 q q q 2 • • • q q q n-1 and its time-shifted counterpart Y Y Y = q q q 2 q q q 3 • • • q q q n , the minimization problem can be recast in matrix form as minimize

A A A Y Y Y -A A AX X X 2 F .

Its solution is simply the least-squares solution

A A A LS = C C C yx C C C -1 xx , where C C C xx = X X XX X X H and C C C yx = Y Y Y X X X H are sample variance-covariance matrices.
Unfortunately, in fluid mechanics we are often in the situation where m n, i.e. the dimension of the state vector q q q is generally much larger than the number of samples we are able to collect. Hence, the estimates C C C xx and C C C yx of the variance-covariance matrices are not converged, implying that C C C xx is likely to be rank-deficient and its inverse ill-posed. As a consequence, the leastsquare solution A A A LS tends to be overly sensitive to minute details of the data used to train the model and to be not statistically representative of the dynamics one tries to describe.

In order to circumvent this issue, it is common practice in regression analysis to regularize the optimization problem by assuming for instance that the unknown matrix A A A is low-rank. If so, it can be factorized as A A A = P P PQ Q Q H , with P P P and Q Q Q ∈ R m×r two rank-r matrices providing bases for the columnspan and the rowspan of A A A, respectively. Without loss of generality, one can additionally assumes that the columns of P P P form an orthonormal set of m-dimensional vectors, i.e. P P P H P P P = I I I. Under these assumptions, the regularized DMD problem then reads minimize

P P P,Q Q Q Y Y Y -P P PQ Q Q H X X X 2 F subject to rank P P P = r, rank Q Q Q = r,
P P P H P P P = I I I.

Although the above optimization problem is non-convex due to the rank constraint, it belongs to a wider class of problems known as Reduced Rank Regression (RRR) whose properties and closed-form solutions have been studied as early as the mid 1970's by Izenman [34] and recently reviewed by De la Torre [35]. It can be shown that the above norm minimization problem is equivalent to the following trace maximization problem

Q Q Q = C C C -1
xx C C C xy P P P Given P P P and Q Q Q, the low-rank DMD operator can easily be factorized in a second step as

A A A = ΨΛΦ H ,
where Ψ, Φ ∈ C m×r are the so-called left and right DMD eigenmodes, respectively, while Λ ∈ C r×r is a diagonal matrix containing the DMD eigenvalues λ i . These eigenvalues may provide valuable insights into the dynamics of the spatio-temporal coherent structures captured by the DMD eigenmodes. A complete derivation of the optimal solution and of the equivalence between the norm-minimization and trace-maximization problems is provided in Appendix A. It must be noted that most of the DMD variants proposed in the literature [12,13,14,26,27,[START_REF] Kutz | Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems[END_REF][START_REF] Dawson | [END_REF]30,31,32] actually compute only suboptimal approximations of the solution to the rank-constrained DMD problem based on various heuristics. To the best of our knowledge, only Héas & Herzet [36] have proposed an alternative derivation of the optimal DMD solution albeit based on more advanced mathematics than the ones used herein. Note finally that, using an information theoretic point of view, it has been shown by Tegmark [37] for linear time-invariant stochastic dynamical systems that using two low-dimensional latent representations (namely P P P and Q Q Q or Ψ and Φ) enables us to better approximate the causal relationship between q q q k and q q q k+1 than when using a single low-dimensional one as done in classical versions of DMD or when using POD to characterize dynamics rather than statistics.

Sparse identification of nonlinear dynamics

Identifying dynamical systems from data has been a central challenge in mathematical physics. The form of the dynamics is typically either constrained via prior knowledge, as in Galerkin projection, or a particular model structure is chosen heuristically and parameters are optimized to match the data. Simultaneous identification of the model structure and parameters is considerably more challenging as there are combinatorially many possible model structures. The sparse identification of nonlinear dynamics (SINDy) approach [18] bypasses the intractable brute force search by leveraging the fact that many systems may be modeled as dx

x x dt = f f f (x x x)
where x x x ∈ R n is the state vector of our system and the unknown f f f : R n → R n is sparse in the space of possible right-hand side functions.

In order to identify f f f (x x x), time-series data is first collected and formed into the data matrix X X X = x x x(t 1 ) x x x(t 2 ) • • • x x x(t m ) T .

A similar matrix Ẋ X X of time derivatives is formed as well. As a second step, a possibly overcomplete library (or dictionary) Θ of candidate nonlinear functions is constructed, e.g.

Θ (X X X) = 1 1 1 X X X P 2 (X X X) • • • P d (X X X) .
Here P d (X X X) denotes a matrix with columns vectors given by all possible time-series of d th degree monomials in the state x x x, e.g.

P 2 (x 1 , x 2 , x 3 ) = x 2 1 x 1 x 2 x 1 x 3 x 2 2 x 2 x 3 x 2 3 .
Note that any basis function may be included in the library Θ, albeit polynomials have proven to work well for fluid problems. The unknown dynamical system may now be represented in terms of the data matrices as

Ẋ X X = Θ(X X X)Ξ.
Each column Ξ k is a vector of coefficients determining the active terms in the equation governing the dynamics of x k (t). A parsimonious model will provide an accurate model fit with as few non-zero terms as possible in Ξ. Identifying this sparsity pattern may be formulated as the following optimization problem minimize

Ξ k card(Ξ k ) subject to Ẋ X X -Θ(X X X)Ξ k 2 2 ≤ σ C C CΞ k = d d d h h h(Ξ k ) ≤ 0,
where card(Ξ k ) = Ξ k 0 is the cardinality (or 0 norm) of Ξ k , i.e. its number of non-zero entries, and the constraint Ẋ X X -Θ(X X X)Ξ k 2 2 ≤ σ quantifies the fidelity of the model with respect to the data. The linear equality constraint C C CΞ k = d d d and the convex inequality constraint h h h(Ξ k ) ≤ 0 may be used to enforce additional constraints on the identified model based on prior physical knowledge such as energy-preserving quadratic non linearity in [38] for instance. The above minimization problem remains however a combinatorially complex optimization problem. Various convex relaxations to this problem have been proposed over the years, most of them replacing the 0 norm with an 1 norm such as in LASSO regression [39]. Due to the extended set of constraints that will be used in this work, the convex relaxation of this problem based on Thikonov regularization and iteratively thresholded least-squares was found to be sufficient.

Results

This section provides a factual presentation of the results. First, results from the dimensionality reduction using dynamic mode decomposition are presented in § 5.1 while § 5.2 focuses on the identification of the low-order model governing the chaotic dynamics of the thermosyphon.

Low-dimensional embedding

The Navier-Stokes equations are integrated forward in time until a statistical steady state is reached. Once reached, a sequence of 10 000 snapshots of the complete state vector of the fluctuation is collected with a sampling period ∆ T = 0.003 diffusive time units. Such a sampling period enables us to properly sample each oscillation of the flow in the clockwise or counter-clockwise direction using approximately 40 snapshots. It must be noted moreover that, as shown in figure 2, the flow exhibits symmetry with respect to the vertical axis. By taking the symmetric of each collected snapshots, this dataset is thus artificially augmented to 20 000 snapshots. This process enables us to enforce the existing spatial symmetry of the high-dimensional system directly in the training dataset itself.

Dedicated DMD analyses are conducted for the velocity and the temperature field separately. Figure 8 depicts the eigenspectrum of the matrix

C C C yx C C C -1
xx C C C xy for both fields. It can be seen that, as for the Lorenz system, the velocity field is well approximated by a rank 1 model capturing close to 98% of the fluctuation's kinetic energy. Similarly, the temperature field is well approximated by a rank 2 model. Note however that in this case the rank 2 model only captures 76% of the fluctuation's thermal energy. Capturing 98% of the fluctuation's thermal energy would require a model of rank 20, thus significantly increasing the complexity of our overall modeling procedure. The spatial distribution of these three DMD modes are shown in figure 9. While the DMD mode associated to the velocity field is almost invariant in the azimuthal direction, the modes associated to the temperature are anti-symmetric with respect to the vertical or horizontal axis. Interestingly, the features captured by these DMD modes (namely the flow rate, the left-right and top-bottom temperature differences) correspond to the measurements considered by Erhard & Müller [22] to characterize experimentally the chaotic dynamics in a similar thermosyphon. They are moreover consistent with the driving modes derived analytically by Yorke et al. [21] using a Bessel-Fourier decomposition for a simplified set of partial differential equations modeling the same problem.

Figure 10 shows the projection of time-series of the high-dimensional state vector onto the span of these DMD modes. As we had expected, the attractor within this low-dimensional space is strikingly similar to the canonical Lorenz attractor, see figure 17 for a visual comparison. Although not shown, it has the same (x, y, z) → (-x, -y, z) equivariance as the canonical Lorenz attractor. These time-series will form the training dataset for the identification of the low-order model in the next section.

System identification

Let us now turn our attention to the identification of a low-order model able to approximate correctly the dynamics of the DMD measurements presented in the previous section. Note

x(t) z(t) 0 1 2 3 4 5
x(t)

0 1 2 3 4 5
y(t) 0 1 2 3 4 5 t z(t)
Fig. 10 Low-dimensional embedding obtained by projecting time-series of the high-dimensional state vector onto the span of the chosen DMD modes. x(t) denotes the amplitude of the DMD mode associated to velocity while y(t) and z(t) are the amplitudes of the first and second DMD modes associated to the temperature field, respectively. Note that, for the sake of simplicity later on, time-series of x(t), y(t) and z(t) have been standardized such that they all have zero-mean and unit-variance.

that the amplitudes x x x(t) of the DMD modes are mean-centered and normalized to unitvariance to ease the identification process. Additionally, we will assume that the unknown system exhibits quadratic nonlinearities (consistent with the Navier-Stokes equations) so that the library Θ(x x x) of candidate functions for the identification is given by

Θ(x, y, z) = 1 x y z x 2 xy xz y 2 yz z 2 .
Under these assumptions, our unknown model takes the general form

ẋ = a 0 + a 1 x + a 2 y + a 3 z + a 4 x 2 + a 5 xy + a 6 xz + a 7 y 2 + a 8 yz + a 9 z 2 , ẏ = b 0 + b 1 x + b 2 y + b 3 z + b 4 x 2 + b 5 xy + b 6 xz + b 7 y 2 + b 8 yz + b 9 z 2 , ż = c 0 + c 1 x + c 2 y + c 3 z + c 4 x 2 + c 5 xy + c 6 xz + c 7 y 2 + c 8 yz + c 9 z 2 .
With no other assumptions about the dynamics, up to thirty coefficients thus need to be identified from our training dataset.

Deriving the constraints

Assuming that the low-order dynamics are governed by a Lorenz-like process, § 3.3 has introduced the various properties the model would have to exhibit. Because our unknown model is linear-in-parameters, these properties can easily be enforced in the identification step by constraining the coefficients of the model.

Equivariant dynamics :

The first property we will discuss is that of equivariance. Analysis of the phase portraits shows that the dynamics of the thermosyphon embedded into the lowdimensional DMD space appear to be equivariant with respect to the change of coordinates

(x, y, z) → (-x, -y, z).
Denoting by

γ γ γ =   -1 0 0 0 -1 0 0 0 1  
the operator performing this change of coordinates, the yet-to-be identified unknown system thus needs to satisfy

γ γ γ ẋ x x = f f f (γ γ γx x x) .
for the dynamics to be equivariant. Given our choice for the library Θ(x, y, z), a dynamical system equivariant under the action of γ γ γ would take the form ẋ = a 0 x + a 1 y + a 3 xz + a 4 yz,

ẏ = b 0 x + b 1 y + b 3 xz + b 4 yz, ż = c 0 + c 1 z + c 2 x 2 + c 3 xy + c 4 y 2 + c 5 z 2 .
Explicitly enforcing the equivariance constraint into the structure of the model thus reduces the number of coefficients to be identified from thirty down to only fourteen.

Dissipative dynamics :

The second property that needs to be verified is that of dissipative dynamics. As shown in § 3.3, this property reads

∇ • f f f (x x x) < 0 ∀x x x.
Starting from the equivariant model, the divergence of the unknown right-hand side function is given by

∇ • f f f (x x x) = a 0 + b 1 + c 1 + (a 3 + b 4 + c 5 ) z.
In order for our model to have dissipative dynamics in the whole phase space, the unknown coefficients must satisfy the following constraints

a 0 + b 1 + c 1 < 0 a 3 + b 4 + 2c 5 = 0.
Note that this set of constraints couples all of the dynamical equations thus forcing us to identify all of them at once rather than sequentially.

Energy-preserving nonlinearities : The quadratic nonlinearities of the incompressible Navier-Stokes equations do not create nor dissipate energy but simply scatters it among the different length scales. As shown in [38] this property can also be enforced in the identification process using linear equality constraints. Starting from the equivariant model, the equation governing the energy E = x 2 +y 2 +z 2 /2 needs to read

Ė = c 0 z + a 0 x 2 + b 1 y 2 + c 1 z 2 + (a 1 + b 0 )xy
in order to ensure that the quadratic nonlinearity of the unknown model is energy-preserving. As a consequence, the coefficients of the equivariant model need to satisfy

b 3 + c 4 = 0, a 2 + c 2 = 0,
a 3 + b 2 + c 3 = 0 c 5 = 0.
Stable manifold of the conducting state : The last property discussed in § 3.3 stated that the z-axis belongs to the stable manifold of the fixed point associated to the pure conducting state. Assuming that x = y = 0, the equation for z simply reads

ż = c 0 + c 1 z.
The corresponding fixed point is given by z =c 0/c 1 . Additionally, this fixed point is stable provided that c 1 < 0.

This linear inequality constraint is the last constraint we will enforce in the identification process.

Identifying the system

Convex relaxation : As discussed in § 4.2, the optimization problem to be solved for the identification is technically a combinatorial problem which becomes rapidly intractable as the number of candidate functions in the library Θ(x x x) increases. In order to bypass this issue, a convex relaxation of the problem based on Thikonov regularization is used in the present work. The corresponding constrained convex minimization problem then reads minimize

Ξ 3 ∑ k=1 Θ(x x x)ξ k -ẋ x x k 2 2 + λ ξ k 2 2 subject to C C CΞ = 0 H H HΞ ≤ 0. (2)
with Ξ = ξ 1 ξ 2 ξ 3 T the vector of unknown coefficients for the three equations and λ a hyper parameter to be optimized. Note that the 2 -regularization used herein differs from the sparsity-promoting approach classically used in SINDy [18,40,41,42,43,44,38,45]. Although Thikonov regularization is not per say a sparsity promoting regularization, the admissible structure of the model is already sufficiently restricted by the physics-derived constraints discussed in the previous section and this simple approach was found to be sufficient.

Model selection : Solving the optimization problem (2) gives rise to a family of candidate models parameterized by λ from which a single one needs to be selected. To do so, each low-order model from this family is integrated forward in time until a statistical steady state is reached and the resulting empirical variance-covariance matrix Σ is compared against the variance-covariance matrix Σ obtained from the training dataset. This comparison relies on the Kullback-Leibler divergence between two symmetric positive-definite matrices. It is defined as

KL Σ Σ = 1 2 Tr Σ-1 Σ -n + ln | Σ| |Σ|
where n is the number of degrees of freedom (here 3) and |Σ| denotes the determinant of said matrix. Figure 11 shows the evolution of this metric as a function of the regularization parameter λ for the two classes of models, namely unconstrained and physics-constrained models. Enforcing prior physical knowledge into the identification procedure leads to orders of magnitude more accurate models. coefficients that need to be identified without prior physical knowledge. Constrained models correspond to models identified using all of the physics constraints discussed in § 5.2.1. Identified model : Based on figure 11, the constrained model obtained for λ 10 -2 is considered the best model. Its coefficients are then refitted using an extended least squares procedure [START_REF] Billings | Nonlinear System Identification: NARMAX Methods in the Time, Frequency, and Spatio-Temporal Domains[END_REF] to remove any left over correlation in the residuals in order to obtain unbiased estimates. Rounding its coefficients to two significant digits, this model reads ẋ = -76.08x + 88.39y ẏ = 20.76x -41.49xz -4.19y ż = 41.49xy -43.67 -17.31z. By design, it verifies all of the physical constraints discussed in § 5.2.1. Except for the constant term in the z-equation, the structure of this model is moreover identical to that of the Lorenz system, see § 3. Figure 12 provides a qualitative comparison of the true strange attractor and the one of the identified model along with representative time series of the various degrees of freedom. Good agreement in the "eye-ball norm" can be observed. More quantitative analyses and discussion about the accuracy of the identified model are postponed to § 6.2.

x(t) z(t) 0 1 2 3 4 5
x(t) DNS Identified model 0 1 2 3 4 5
y(t) 0 1 2 3 4 5 t z(t)

Discussion

After having presented the results in § 5, this section provides a more in-depth discussion about the model's properties as well as its interpretation in terms of physical mechanisms. Discussion pertaining to dynamic mode decomposition is presented in § 6.1 while properties of the identified low-order model are discussed in § 6.2.

On the dynamic mode decomposition

In this work, choice has been made to learn a low-dimensional embedding for the velocity and temperature fields separately. Based on the eigenvalue distribution of symmetric positive definite matrices, a rank 1 model was found to be sufficient for the velocity field while a rank 2 model (albeit less accurate) was deemed enough for the temperature field. Following [46], one could have alternatively defined the inner product inducing the total energy as q q q i , q q q j = 1

2 Ω u u u i • u u u j + γθ θ θ i • θ θ θ j dΩ (3)
and learn the embedding for the complete state vector at once. Here, γ weights the kinetic energy and the thermal one. Depending on its value, different embeddings could be obtained. When γ is set to Ra Pr and i = j, the first term describes the kinetic energy of the fluctuation while the second one corresponds to its non-dimensional potential energy. Figure 13 (a) shows the eigenvalue distribution of the symmetric positive definite matrix C C C y y yx x x C C C -1 x x xx x x C C C x x xy y y when the whole state vector is considered at once for γ = Ra Pr while figure 13 (b) shows how the total energy is distributed between kinetic and potential energy for the leading eigenvectors of this symmetric positive definite matrix. Looking at these two figures, one can observe that a rank 1 model captures close to 95% of the linearly predictable evolution of the fluctuation's total energy. Additionally, almost all of this total energy consists in kinetic energy. On the other hand, the total energy of the second leading eigenvector mainly consists in potential energy. This clear distinction between the distribution of energy for the first two eigenvectors supports our decision to learn a separate embedding for the velocity and the temperature.

Figure 14 (a) depicts the azimuthal velocity component and temperature field of the DMD modes when a rank 3 model is considered while figures 14 (b) and (c) provides a comparison of the low-dimensional representation obtained from learning jointly or separately the embedding of the state vector, respectively. Because they are highly linearly correlated, the real part of the first DMD mode includes both the leading coherent structure associated to the velocity and the left-right temperature fluctuation. As a consequence, these two structures are fused into a single degree of freedom when the embedding for the velocity and the temperature is learned jointly. With the exception of its velocity component, the third DMD mode captures as before the leading coherent structure associated to the vertical temperature fluctuation. Figures 14(b) and (c) provide a comparison of the resulting lowdimensional representations. While the projection onto the (x, z) plane appears similar for both embeddings, the projection onto the other planes are different. Despite our efforts, we have not been able to successfully identify a low-order model when learning the embedding for the velocity and temperature jointly. This failed attempt highlights the importance of choosing judiciously the embedding when aiming for the identification of an interpretable low-order model and seems to suggest that learning a low-dimensional embedding for each physical quantity separately might provide more useful information into the dynamics of the system than when learning the embedding for the whole state vector at once. This might be better understood by realizing that, given the state vector q q q = v v v θ θ θ

T , the DMD operator may be partitionned as

v v v k+1 θ θ θ k+1 = A A A v v vv v v A A A v v vθ θ θ A A A θ θ θ v v v A A A θ θ θ θ θ θ v v v k θ θ θ k .
Performing a component-based DMD analysis thus enables us to obtain different (and possibly more accurate) low-rank approximations of each block while also providing a finer-grain description of the interplay between velocity and temperature. To date, this observation however lack theoretical justifications.

6.2 On the identified low-order model

How accurate is it?

While figure 12 shows good qualitative agreements between the dynamics of the thermosyphon embedded within the low-dimensional DMD space and that of the identified model, let us now try to characterize this agreement more quantitatively. Figure 15 provides a comparison of the probability distribution functions obtained from direct numerical simulation of the Navier-Stokes equations and those obtained from the identified low-order (a) DMD modes obtained using the total energy inducing inner product

x(t) y(t) x(t) z(t) y(t) z(t) (b)
Embedding for the velocity and temperature learned jointly

x(t) y(t) x(t) z(t) y(t) z(t)
(c) Embedding for the velocity and temperature learned separately 16 Comparison of the autocorrelation function and time-scale distributions obtained from direct numerical simulation (light gray) and the low-order model (black). For the direct numerical simulation, these statistical quantities have been computed using the time-series of the flow rate while that of x(t) (associated to the velocity DMD mode) has been used for the low-order model. Both autocorrelation functions have been computed for the corresponding telegraph signal. model. Given that the low-order model only has three degrees of freedom, excellent agreement is obtained for all three probability density functions. A small discrepancy is however observed for small values of x(t) and y(t). It corresponds to situations where the flow is about to switch from clockwise to anti-clockwise direction (or vice-versa). During these events, the amplitude of the DMD modes associated to the velocity and the left-right temperature difference are close to zero. Hence, a plausible explanation for this mismatch is that, during these events, the dynamics of the flow are dominated by other coherent structures not accounted for by the present model. The overall dynamics of the system are nonetheless correctly captured as assessed by figure 16. Using the same symbolic dynamics approach as in § 2.3, it compares (a) the autocorrelation function of the telegraph signal and (b) the time-scale distribution obtained from direct numerical simulation of the Navier-Stokes equations and by our low-order model. Remarkable agreement is obtained for the autocorrelation function. Reasonable agreement is also obtained regarding the time-scale distribution. For all intent and purposes the low-order model identified in the present work thus captures remarkably accurately the statistical and dynamical properties of the flow.

Interpretability and physical processes captured by the model

Over the past few years, interpretability has been an overarching goal in machine learning, in particular when applied to modeling of physical systems. Because it relies on the joint use of dynamic mode decomposition and SINDy, every single mode of the low-dimensional embedding and term of the low-order model identified in the present work can be understood as modeling a given coherent structure or physical mechanism at play in the thermosyphon. Additionally, a sequence of linear change of variables and re-scaling time enables us to recast the identified model into the Lorenz system ẋ = σ (yx) (18.14, 34.75, 4.13), thus effectively reducing the number of parameters from seven down to only three and facilitating their interpretability and analysis.

ẏ = ρx -y -xz ż = xy -β z with (σ , ρ, β ) =
Analysis of the above model shows that it has three fixed points. The trivial one, given by (x, y, z) = (0, 0, 0), corresponds to the pure conducting state (i.e. only temperature diffusion and no net flow through a cross section of the thermosyphon). As expected, this particular fixed point is moreover linearly unstable with its unstable eigenvector corresponding to the linear seed of a convection cell. The other two fixed points are symmetric of one another and correspond to the clockwise and anti-clockwise configurations. Interestingly, given the set of identified parameters, these fixed points appear to be very slightly stable, thus suggesting that the dynamics observed in the thermosyphon do not correspond to full-blown chaos but only transient chaos. It is unclear at the present time whether this is an actual feature of the flow or only a limitation of the identification procedure used in this work. It must be noted however that simulations of the identified model show that this transient chaos can nonetheless sustains itself for a period of time larger than hundred times the complete duration of the direct numerical simulation used to generate our training data. Moreover, increasing the training data with another set of 10 000 snapshots and re-running the whole identification procedure leaves our results largely unchanged. In the absence of further evidence, one cannot conclude about the nature of the dynamics (chaos or only transient chaos) of the Navier-Stokes equations for the set of physical parameters considered herein.

Because it reduces to the canonical Lorenz system, the various terms in the identified can be interpreted as modeling similar physical mechanisms as those observed in the canonical Rayleigh-Bénard convection. For instance, in the governing equation for x(t), the term -σ x serves to model momentum diffusion due viscous dissipation while σ y captures the influence of buoyancy on the flow. Similarly for the equation governing y(t), the term -y captures the diffusion of left-right temperature fluctuation while the terms ρx and -xz serve to model the advection of the mean temperature profile and vertical temperature fluctuation, respectively. Finally, the terms xy and -β z in the last equation model the advection of the left-right temperature fluctuation and the diffusion of vertical temperature fluctuation. Regarding the parameters of the model, it is well known from the Lorenz system that ρ is the ratio of the Rayleigh number Ra over the critical Rayleigh number Ra c corresponding to the onset of the convection cells. Similarly, the parameters σ and β can be directly related to the Prandtl number and the aspect ratio of the convection cells. Given that no prior knowledge about the system other than the physical constraints it needs to satisfy had been included in the identification process, the value ρ = 34.75 is in excellent agreement with the expected ratio of Ra over Ra c . Indeed, preliminary computations of the complete bifurcation diagram of the high-dimensional system shows that the critical Rayleigh number for the onset of the steady convection cells is Ra c = 495. For the Rayleigh number considered in this work, this would thus correspond to ρ = 34.34, i.e. only 1% different than the value identified. It is thus quite remarkable that a model identified at Ra = 17 000 is able to accurately predict a fundamental property of the high-dimensional such as a bifurcation happening at a fairly different value of the control parameter. Note however that, since a single operating condition has been considered herein, this excellent agreement might be only coincidental and further analyses are required to confirm it.

Conclusion and perspectives

Identifying accurate yet interpretable low-order model from data has been an overarching goal in mathematical physics which has gained a renewed interest over the past decade with the advances made in dimensionality reduction and sparsity-promoting regression techniques. In this work, we have illustrated how such techniques, namely dynamic mode decomposition [13] and SINDy [18], could be used to identify a low-order model of the chaotic thermal convection taking place in an annular thermosyphon at moderate Rayleigh numbers. Particular emphasis has been put on the necessity for the low-dimensional embedding and the identified system to comply with prior physical knowledge. To do so, the construction of the low-order model proposed herein closely follows the derivation of the canonical Lorenz system from the Oberbeck-Boussinesq equations modeling the Rayleigh-Bénard convection between two infinite flat plates. Given the simplicity of the identified model compared to the numerical simulation of the original high-dimensional system, excellent qualitative and quantitative agreements have been obtained as assessed by the similar structure of the strange attractor onto which the dynamics evolve in both situations as well as various other statistical quantities.

From a physical point of view, this work further confirms a strand of evidence [21,22,23] that the chaotic dynamics resulting from unstable temperature stratification in annular thermosyphons may be correctly captured using a low-order model as simple as the canonical Lorenz system. Numerous questions however remain unanswered, the most important ones being -Are the conclusions drawn from this study applicable to other configurations? Is the Lorenz model a universal low-order model for other unstable temperature stratifications (e.g. Rayleigh-Bénard flows in rectangular enclosures [47,48]) at low or moderate Rayleigh numbers? -To what extent does the Lorenz system, combined with a DMD embedding, provides a good approximation of the dynamics as the Rayleigh number increases? How should one modify the model in order to account for turbulence once it sets in?

While it may be relatively easy to provide an answer to the first questions by running the same analyses as those conducted herein for other geometries, the answer to the second questions is highly dependent on how compressible (from a dimensionality point of view) the dynamics of the high-dimensional system are once turbulence sets in. Finally, it should be emphasized that our ability to identify low-order models from data stemming from a high-dimensional system with various practical applications opens new possibilities. In particular, both dynamic mode decomposition and SINDy can easily be extended to include external input [49,41], enabling us to identify accurate low-order models for control purposes. These models could moreover be used in conjunction with optimal sensor placement [50,51] to enable real-time monitoring and estimation of the system. Using tool from model predictive control, such low-order models could be used to strengthen or decrease the chaotic dynamics depending on the applications. A first step toward this goal has been achieved by Kaiser et al. [44].

A Derivation of the optimal solution to the DMD problem

A.1 Equivalence between norm-minimization and trace-maximization formulation of the DMD problem

This section aims at deriving the equivalence between the norm-minimization and trace-maximization formulations of the DMD problem and the corresponding solution. Given two data sequences X X X and Y Y Y related by y y y k = f f f (x x x k ), the aim of DMD is to find the low-rank operator A A A that best approximates the possibly nonlinear function f f f in the least-squares sense. As discussed in § 4.1, this can be formulated as the following rank-constrained minimization problem minimize

P P P,Q Q Q Y Y Y -P P PQ Q Q H X X X 2 F subject to rank P P P = rank Q Q Q = r P P P H P P P = I I I, (4)
where we use the fact the unknown low-rank operator can be factorized as A A A = P P PQ Q Q H . Introducing the cost function

J (P P P, Q Q Q) = Y Y Y -P P PQ Q Q H X X X 2 F = Tr (Y Y Y -P P PQ Q Q H X X X)(Y Y Y -P P PQ Q Q H X X X) H
for the sake of notational simplicity, the solution (P P P * , Q Q Q * ) of Eq. ( 4) is given by the stationary points of J (P P P, Q Q Q), i.e. ∂ J

∂ P P P = -2Y Y Y X X X H Q Q Q + 2P P PQ Q Q H X X XX X X H Q Q Q = 0 ∂ J ∂ Q Q Q = -2X
X XY Y Y H P P P + 2X X XX X X H Q Q QP P P H P P P = 0

Introducing the orthogonality constraint on P P P in the second equation above yields the following expression for

Q Q Q Q Q Q H = P P P H C C C yx C C C -1 xx , where C C C yx = Y Y Y X X X H and C C C xx = X X XX X X
From this point, it is straightforward to show that the above problem can be rewritten as

Finally, recalling that the trace of a matrix is the sum of its eigenvalues, basic linear algebra is sufficient to prove that the solution of the above optimizations problem is given by the first r leading eigenvectors of the symmetric positive definite matrix C C C yx C C C -1 xx C C C xy . Once P P P has been computed, the matrix

Q Q Q is simply given by Q Q Q H = P P P H C C C yx C C C -1
xx as stated previously. Interestingly, it can be noted that if X X X is a full rank m × m matrix (i.e. as many linearly independent samples as degrees of freedom), then

C C C yx C C C -1 xx C C C xy = Y Y Y X X X H X X XX X X H -1 X X XY Y Y H = Y Y YY Y Y H = C C C yy .
In this particular case (hardly encountered in fluid dynamics due to the memory footprint of storing m snapshots of a m-dimensional vector), the columns of P P P are simply the first r leading eigenvectors of the variancecovariance matrix C C C yy , i.e. the POD modes of the output matrix Y Y Y . Alternatively, if X X X is a skinny m × n matrix (with m > n) or a rank-deficient m × m matrix, one can introduce its economy-size singular value decomposition X X X = U U U X X X Σ X X X V V V H X X X . The columns of P P P are then given by the first r left singular vectors of the matrix Y

Y YV V V X X X V V V H X X X where V V V X X X V V V H X X X
is the projector onto the rowspan of X X X.

A.2 Computing the DMD modes and eigenvalues from the low-rank factorization of A A A

It must be noted that the columns of P P P are not the so-called DMD modes Ψ although they span the same subspace. The DMD modes can nonetheless be easily computed from the low-rank factorization A A A = P P PQ Q Q H of the DMD operator. The i th DMD mode is solution to the eigenvalue problem

Ψ i λ i = A A AΨ i .
Since Ψ i needs to be in the columnspan of A A A, it can be expressed as a linear combination of the columns of P P P, i.e. Although A A A is a high-dimensional m × m matrix, Q Q Q H P P P is a small r × r matrix whose eigenvectors b b b i and eigenvalues λ i can easily be computed using direct solvers thus enabling the computations of the eigenvalues and eigenvectors of A A A at a reduced cost. Similarly, the adjoint DMD modes Φ i are solution to the adjoint DMD eigenvalue problem

Φ i λi = A A A H Φ i ,
where the overbar denotes the complex conjugate. These now live in the rowspan of A A A and can thus be expressed as

Φ i = Q Q Qc c c i
where once again c c c i is a small r-dimensional vector. The corresponding low-dimensional eigenvalue problem then reads c c c i λi = P P P H Q Q Qc c c i with P P P H Q Q Q a small r × r matrix. It must be noted however that, from a practical point of view, converging the adjoint DMD modes Φ i may require significantly more data than needed to the converge the direct DMD modes Ψ i .

B Broomhead-King embedding for attractor reconstruction

Attractor reconstruction from a single time-series has been instrumental in highlighting the connection between the dynamics of the thermosyphon and that of the canonical Lorenz system (see figure 3 and figure 4). The aim of this section is to briefly introduce inexperienced readers to the Broomhead-King attractor reconstruction technique used herein. For more details, please refer to the original paper [20]. For the sake of simplicity and reproducibility, we will apply this attractor reconstruction technique using the Lorenz system ẋ = σ (yx) ẏ = x (ρz)y ż = xyβ z with (σ , ρ, β ) = (10, 28, 8 /3) so that it exhibits chaotic dynamics. The properties of this dynamical system have been discussed in § 3. Figure 17 depicts its well-known strange attractor along with representative timeseries of x(t), y(t) and z(t). In the rest of this appendix, we will work exclusively with the time-series of x(t).

The whole dataset consists of the evolution of x(t) from t = 0 to t = 400 with a sampling period ∆t = 0.001 resulting in 400 000 samples. Although the Takens embedding theorem [19] provides theoretical guarantees to reconstruct a multidimensional attractor from a single time-series, this approach might break when applied to real data that may be contaminated by measurement noise or it might be quite sensitive on the sampling period used to acquire the signal. Additionally, once the lag τ has been determined, the number of time-lagged versions of x(t) needed to reconstruct the attractor needs to be determined using techniques based nearest neighbors [52,53,54] which might provide conflicting estimates of the attractor's dimension and different reconstructions. As to overcome these limitations, Broomhead & King [20] have proposed instead to reconstruct the attractor based on the singular value decomposition of the Hankel matrix constructed from time-lagged version of x(t)

H H H = 1 √ N      x 1 x 2 x 3 • • • x d x 2 x 3 x 4 • • • x d+1 . . . . . . . . . • • • . . . x N-d+1 x N-d+2 x N-d+3 • • • x N     
where x k = x(t k ), N is the total length of the time-series (i.e. 400 000 in our case) and d is the maximum number of time-lagged considered. The choice of d will be discussed later on. Decomposing this Hankel matrix as H H H = U U UΣV V V H , the dimension of the attractor can then be inferred from the distribution of the singular values σ i . The columns of U U U then provide an orthonormal basis for the embedding space while the columns of V V V describe the evolution of the system within this embedding space, thus enabling the reconstruction we aimed for. The choice of the number of lags d needed to construct the Hankel matrix is of crucial importance as it will determine the window length τ w = d∆t used to embed the dynamics which, in turn, will determine the structure of the singular value spectrum. In practice, the window length τ w needs to be large enough to average out the possible noise contaminating the data. Note however that, in the limit τ w → ∞, the singular value decomposition of the Hankel matrix converges to the discrete Fourier transform of x(t). Chaotic systems being characterized by a continuous spectrum, an infinite number of singular components would thus be needed to reconstruct the signal, grossly overestimating the dimension of the attractor. In their paper, Broomhead & King [20] recommend to use a window length τ w of approximately one tenth of the period of oscillation about the unstable saddle-foci. In this work, this window length has been chosen as one twentieth of this oscillating period, even though choosing one tenth does not change qualitatively (nor quantitatively) the results presented. Note that the same rationale has been used when applying this technique to the data from the thermosyphon. For our choice of parameters, figure 18 provides the temporal evolution of the system within the embedding space and the corresponding reconstruction of the attractor. The theoretical properties of such time-delay embedding of attractors have been recently studied by Kamb et al. [55].

Fig. 1

 1 Fig.1Left: Sketch of the geometry considered. The temperature of the upper walls of the thermosyphon (dashed lines) is set to θ = 0 while it is set to θ = 1 on the lower walls (solid lines). The gravity is acting along -e e e y . Right: Definition of the various non-dimensional parameters defining the problem. Throughout this manuscript, the Rayleigh number is set to 17 000 while the Prandtl one is set to 5.

Fig. 2

 2 Fig. 2 Temperature, radial velocity and azimuthal velocity fields at various instants of time. The same colorscale has been used for each row.

  Fig. 3 (a) Time-series of the flow rate extracted from a direct numerical simulation of the annular thermosyphon at Ra = 17 000 and Pr = 5. The whole time-series extends from t = 0 to t = 650. Only a subset is shown. (b) Empirical probability distribution of the flow rate m(t). (c) Power spectral density of the timederivative ṁ(t) of the flow rate. (d) Broomhead-King reconstruction of the underlying strange attractor from the time-series shown in (a). The temporal window used for this reconstruction is set to τ = 0.008 nondimensional time-units corresponding to 32 measurements of the flow rate. Each panel of figure (b) shows the projection of the strange attractor onto the planes (v 1 , v 2 ), (v 1 , v 3 ) and (v 2 , v 3 ), respectively, where v i is the i th Broomhead-King coordinate.

  Fig. 5 (a) Mean flow computed over 650 non-dimensional time units. The light gray line in the temperature plot depicts the θ (x, y) = 0.5 isocontour. For the azimuthal velocity field, it depicts the u a (x, y) = 0 isocontour. (b) Diagonal entries of the Reynolds stress tensor corresponding to (d) θ θ , (e) u r u r and (f) u a u a where prime variables denote fluctuating quantities.

Fig. 6 Fig. 7

 67 Fig. 6 Conditionally averaged mean flow ū u u CCW (x x x) computed over 650 non-dimensional time units. The light gray line in the temperature plot depicts the θ (x, y) = 0.5 isocontour. For the azimuthal velocity field, it depicts the u a (x, y) = 0 isocontour.

  maximize P P P Tr P P P H C C C yx C C C -1 xx C C C xy P P P subject to rank P P P = r, P P P H P P P = I I I. The P P P matrix is thus formed by the first r eigenvectors of the symmetric positive-definite matrix C C C yx C C C -1 xx C C C xy . Note that, if and only if X X X is a full-rank m × m matrix, these reduce to the first r left singular vectors of the output matrix Y Y Y (i.e. the POD modes of Y Y Y ). Once the P P P matrix has been identified, the low-rank Q Q Q matrix is given by

Fig. 8

 8 Fig. 8 Eigenspectrum of the matrix C C C yx C C C -1 xx C C C xy for the velocity (blue circles) and the temperature (gray squares). For the velocity field, a rank 1 DMD model captures 98% of the fluctuation's kinetic energy. For the temperature field, a rank 2 DMD mode captures 76% of the fluctuation's thermal energy. For the latter, a DMD model of rank 20 would be needed to capture 99% of the thermal energy.

Fig. 9

 9 Fig. 9 Spatial distribution of the DMD modes retained for our reduced-order model. Projection of the dynamics onto the span of these DMD modes is shown in figure 10.

Fig. 11

 11 Fig.11Evolution of the Kullback-Leibler divergence of the variance-covariance matrices as a function of the regularization parameter λ . Unconstrained models correspond to the most general models with up to 30 coefficients that need to be identified without prior physical knowledge. Constrained models correspond to models identified using all of the physics constraints discussed in § 5.2.1.

Fig. 12

 12 Fig. 12 Comparison of the dynamics of the chaotic thermosyphon embedded in the low-dimensional DMD space (light gray) and those of the identified low-order model (black). Left panel shows a projection of the Lorenz-like strange attractor onto the (x, z) plane. Right panel shows typical time series of the different degrees of freedom.

Fig. 13 (

 13 Fig. 13 (a) Eigenspectrum of the matrix C C C yx C C C -1 xx C C C xy from the DMD optimization problem when the embedding for the velocity and the temperature is learned jointly using the inner product (3). (b) Energy budget (i.e. fraction of the kinetic and thermal energy forming the total energy) for the leading eigenvectors of the matrix C C C yx C C C -1 xx C C C xy . Note that the DMD modes considered in this section are a linear combination of the first three eigenvectors.

Fig. 14 (Fig. 15

 1415 Fig.14(a) Azimuthal velocity (top row) and temperature (bottom row) fields of the DMD modes associated to the rank 3 model obtained using the inner product inducing the total energy. The first two columns correspond to the real and imaginary part of the complex-conjugate DMD mode while the last column corresponds to the real-valued DMD mode. To be compared with figure9. (b) Low-dimensional representation obtained after projecting the training dataset onto the span of the modes above. (c) Low-dimensional representation used in this work obtained after learning the embedding for the velocity and the temperature separately.

2 F

 2 H are the empirical variance-covariance matrices introduced in § 4.1. Inserting this expression for Q Q Q in Eq. (4) yields minimizeP P P Y Y Y -P P PP P P H C C C yx C C C -1xx X X X subject to rank P P P = r P P P H P P P = I I I,

minimize

  Tr (C C C yy ) -2Tr P P P H C C C yx C C C -1xx C C C xy P P P subject to rank P P P = r P P P H P P P = I I I.

( 6 )

 6 Tr (C C C yy ) being constant, we thus established the equivalence between the norm minimization problem (4) and the following trace maximization problem maximize P P P Tr P P P H C C C yx C C C -1 xx C C C xy P P P subject to rank P P P = r P P P H P P P = I I I.

  Ψ i = P P Pb b b i where b b b i is a small r-dimensional vector. Introducing this expression into the DMD eigenvalue problem yields after some simplifications b b b i λ i = Q Q Q H P P Pb b b i .

Fig. 17

 17 Fig.17Projection of the Lorenz attractor on the (x, z) plane (left) and corresponding time-series of x(t), y(t) and z(t) (right). The parameters of the system are chosen as (σ , ρ, β ) = (10, 28, 8 /3).

Fig. 18

 18 Fig.18Projection of the Lorenz attractor reconstructed using the Broomhead-King technique on the (v 1 , v 2 ) plane (left) and corresponding time-series of v 1 (t), v 2 (t) and v 3 (t) (right). A window length τ w = 0.032, corresponding to roughly one twentieth of the oscillating period around the unstable saddle-foci, has been chosen.

A video of these dynamics is available online at https://tinyurl.com/y55buvnc

Acknowledgements This work is adapted from a contribution to the long program Machine Learning for Physics and the Physics of Learning, organized by the Institute for Pure and Applied Mathematics (IPAM) at UCLA (Los Angeles, USA) in 2019 and made possible thanks to the financial support of IPAM, which is supported by the National Science Foundation (NSF Grant No. DMS-1440415). The author would also like to gratefully thank Steven Brunton, Jared Callahan, Kathleen Champion, Onofrio Semeraro, Alessandro Bucci and many others for the fruitful dicussions on this topic.


Bixin migration kinetic to ethanol 95% v/v (used as a fat food simulant) was evaluated over 65 days at 40 °C. Quantification was made by spectroscopy and the overall mass transfer coefficient (K) of bixin was obtained using the following equation:

The effect of the thermal treatments over films color stability was evaluated using a Konica Minolta colorimeter through L*(luminosity), a* (red-green), and b* (yellow-blue) parameters. The total color differences (∆E) caused by thermal processes when compared to casting films was calculated using the following equation

∆𝐸 = (𝐿 * -𝐿 * !"#$ ) % + (𝑎 * -𝑎 * !"#$ ) % + (𝑏 * -𝑏 * !"#$ ) % where L*,

Results and Discussion

Bixin was extracted from annatto seeds (Bixa orellana) with dichloromethane. Polylactide (NaturePlast PLE 001) was used to develop the active films, where two formulations (F-without plasticizer; A-with plasticizer) were produced by different methods. Formulations and methods are shown below. The bixin loss due to thermal processing of films was evaluated by films dissolution in dichloromethane, followed by bixin quantification by spectroscopy (λ=467 nm; A 1%= 2826). ü The presence of ATBC increased in 10% the bixin loss in TP films.

Bixin (%) ATBC (%) Film

Considering PLA is frequently used on the encapsulation of diversal bioactive compounds -presenting a protective effect -its replacement by ATBC might have reduced bixin stability in films matrix.

ü The presence of ATBC greatly increased bixin migration rate from casting and TP films, which has been ascribed to the increased polymer chain mobility caused by the plasticizer, favoring the active compound release.

ü IMTP films presented the highest overall mass transfer coefficient (K) of bixin to the food simulant.

The cleavage of PLA long chains caused by internal mixer shear forces might facilitate bixin diffusion through the polymeric matrix, increasing its release rate to the food simulant.

The presence of ATBC and mixing shear forces affected bixin stability and migration rate. Since antioxidants' concentration and release rate greatly affect the performance of antioxidant active packaging, these findings are of great importance and are currently being explored on the extension of fatty food products shelf-life. The incorporation of antioxidants into a polymeric structure to produce antioxidant active materials is an alternative way to extend the shelflife of fatty food products. Even tough polylactide (PLA) present satisfactory mechanical properties, it requires improvements to be used as an antioxidant material for food packaging.

Bixin, a carotenoid with hydrophobic character, presents coloring and antioxidant properties and hence is a potential antioxidant additive for PLA. Considering bixin thermal low stability, films produced with bixin and the plasticizer acetyl tributyl citrate (ATBC) were produced by different thermal processes. The effect of thermal treatments over films color changes and bixin stability was evaluated, also as its effect on bixin's kinetic of migration.

  a* and b* were measured in TP and IMTP films. POLYLACTIDE FILMS PRODUCED WITH BIXIN AND ATBC: EFFECT OF THERMAL TREATMENT ON BIXIN STABILITY AND MIGRATION KINETIC Contact Sandra Domenek: sandra.domenek@agroparistech.fr Liana Stoll: lianastoll@hotmail.com JRU 1145 Food Processing and Engineering between INRA and AgroParisTech, Group "Interactions between Materials and Media in Contact", 1 avenue des Olympiades, 91300 Massy, France This study is part of a collaborative research project between AgroParisTech and the Federal University of Rio Grande do Sul (UFRGS), funded by the brazilian program PDSE -CAPES 2019

  Stoll et al. Poly(acid lactic) films with carotenoids extracts: Release study and effect on sunflower oil preservation. Food Chem 281 (2019) 213-221 Different letters in the same column indicate significant differences (p < 0.05) among the treatments

Table 1 .

 1 PLA films formulations and method of obtention(Stoll et al., 2019)

	(1)

Table 2 .

 2 Color parameters, bixin loss and bixin overal migration rate (K) of PLA films Figure 1. Bixin migration from films to ethanol 95% (v/v) Liana Stoll, Alessandro de Oliveira Rios; Sandra Domenek