An Efficient Data-Driven Multiscale Stochastic Reduced Order Modeling Framework for Complex Systems

Many complex turbulent nonlinear systems in nature have the following abstract mathematical structure [6, 1, 9, 47],

where $\mathbf{u}\in\mathbb{C}^{N}$ is the multidimensional state variable and $t$ is time. In (1), $L\mathbf{u}$ and $D\mathbf{u}$ represent linear dispersion and dissipation effects, respectively, where $L$ is skew-symmetric and $D$ is negative-definite. The nonlinear effect is introduced through an energy-conserving quadratic form with $\mathbf{u}\cdot B\left(\mathbf{u},\mathbf{u}\right)=0$. The system is also subject to external forcing effects that contain a deterministic component, $\mathbf{F}\left(t\right)$, and possibly a stochastic disturbance as well, which is assumed to be a Gaussian random noise, $\bm{\sigma}\left(t\right)\dot{\mathbf{W}}\left(t\right)$, for simplicity. The model in (1) can be understood as discretizing a (stochastic) partial differential equation by projecting it onto a given set of basis functions $\{\bm{\varphi}_{1},\bm{\varphi}_{2},\ldots,\bm{\varphi}_{N}\}$, resulting in a high-dimensional system, where the stochastic forcing accounts for the truncation error. Due to the intrinsic turbulent behavior, the system in (1) is highly nonlinear and the state variable $\mathbf{u}$ can possess strongly non-Gaussian statistics in both the marginal and joint PDFs.

Decompose the state variable $\mathbf{u}$ into $\mathbf{u}=(\mathbf{v}^{\mathtt{T}},\mathbf{w}^{\mathtt{T}},\mathbf{z}^{\mathtt{T}})^{\mathtt{T}}$ with $\cdot^{\mathtt{T}}$ being the vector transpose, where $\mathbf{v}\in\mathbb{C}^{N_{1}}$, $\mathbf{w}\in\mathbb{C}^{N_{2}}$, $\mathbf{z}\in\mathbb{C}^{N_{3}}$, and $N_{1}+N_{2}+N_{3}=N$. A typical criterion for such a decomposition is based on the rank of the scales for the state variables. In other words, $\mathbf{v}$ includes the leading few modes corresponding to the large-scale components of the state variables, $\mathbf{w}$ contains the next a few modes describing the medium-scale features, while $\mathbf{z}$ involves the remaining modes that are treated as the small-scale variables. Note that different from many traditional ROM development, where the state variables are decomposed into only the large- and small-scale components, a three-scale decomposition is utilized here. Such a new decomposition facilitates a smoother transition from large- to small-scale variables by means of a set of intermediate-scale variables. It also allows the feedback from medium- to large-scale variables to be explicitly incorporated into the ROM, which is often crucial in enhancing the model accuracy in characterizing energy transfer and intermittency. Finally, the choice of $N_{1},N_{2}$, and $N_{3}$ is case dependent but in a typical scenario $\mathbf{z}$ accounts for the majority number of the modes.

A natural way to build ROMs with a reduced dimension of the coupled system in (1) is to project $\mathbf{u}$ onto its large- and medium-scale components $\mathbf{v}$ and $\mathbf{w}$,

where $\cdot^{*}$ is the conjugate transpose. It is seen that the linear and the quadratic nonlinear terms in (1) have been explicitly expressed in (2). Therefore, $\mathbf{A}^{(\cdot)}_{\cdot}$ and $\tilde{\mathbf{B}}^{(\cdot)}_{\cdot}$ in (2) are all constant matrices and tensors while the two residual terms include the contribution from the variable $\mathbf{z}$ that has been truncated here. The notation in (2) is slightly abused in the sense that the quadratic nonlinear terms between $\mathbf{v}$ and $\mathbf{w}$ are simply written as $\mathbf{v}^{*}\tilde{\mathbf{B}}^{(v)}_{vw}\mathbf{w}$ and $\mathbf{v}^{*}\tilde{\mathbf{B}}^{(w)}_{vw}\mathbf{w}$ while the other parts of the associated nonlinear terms, say $\mathbf{w}^{*}\tilde{\tilde{\mathbf{B}}}^{(v)}_{vw}\mathbf{v}$ and $\mathbf{w}^{*}\tilde{\tilde{\mathbf{B}}}^{(w)}_{vw}\mathbf{v}$, are omitted in (2) for notation simplicity. Note that if the state variables are all real-valued, then $\mathbf{v}^{*}\tilde{\mathbf{B}}^{(v)}_{vw}\mathbf{w}$ and $\mathbf{v}^{*}\tilde{\mathbf{B}}^{(w)}_{vw}\mathbf{w}$ are sufficient to characterize the quadratic nonlinear interactions between $\mathbf{v}$ and $\mathbf{w}$. Since each complex-valued variable can always be rewritten as two real-valued variables, representing its real and imaginary parts, all the state variables in the following will be regarded as real-valued. Similarly, the notation $\cdot^{*}$ hereafter simply denotes the transpose of a real-valued vector. In the development of ROMs, various methods have been proposed to deal with the residual terms. If these residual terms are simply ignored, then the resulting ROM becomes the so-called Garlerkin ROM (G-ROM) [26, 27]. However, it has been shown that the G-ROM often introduces inaccuracy in reproducing even the large-scale dynamics of the FOM. Such a finding implies the importance of building approximate functions to compensate the contribution from the residuals. If these functions depend only on the state variables of the ROM, then they are known as the closures since they allow the ROM to be closed [30, 31, 21, 23, 33]. In most of the existing methods, the closure terms are assumed to be a linear or quadratic function of the entire ROM state variables, which are $(\mathbf{v},\mathbf{w})$ in the setup of (2).

The new ROM to be developed here has several fundamental differences compared with the existing ROMs. First, only the first three terms on the right hand side of (2) are retained in the new ROM framework, which represent the basic dynamics. In other words, the self nonlinear interactions of $\mathbf{w}$ are omitted. One motivation for such a simplification is that $\mathbf{w}$ is a smaller scale variable compared with $\mathbf{v}$ and therefore it is often a faster variable and contains less energy. As a consequence, $\mathbf{w}^{*}\tilde{\mathbf{B}}^{(v)}_{vw}\mathbf{w}$ and $\mathbf{w}^{*}\tilde{\mathbf{B}}^{(w)}_{vw}\mathbf{w}$ are expected to have a weaker role than the other nonlinear terms in (2). Another important reason is related to the development of an efficient data assimilation scheme for the ROM, which will be discussed in Section 2.2. Next, a hybrid closure configuration is introduced to approximate the residual part as well as the contribution of the self nonlinear interactions of $\mathbf{w}$. Here, the closure includes both a set of deterministic functions and a stochastic component. The building block of the deterministic part of the closure in this new ROM framework is very distinct fromthat in the existing methods, which is either linear or quadratic nonlinear in terms of both $\mathbf{v}$ and $\mathbf{w}$. Here, a special type of nonlinear functions are adopted, which are fully nonlinear with respect to $\mathbf{v}$ but are conditional linear with respect to $\mathbf{w}$. The motivation is again that, due to the relatively faster nature of $\mathbf{w}$, a conditional linear function and stochastic terms can effectively approximate the statistical behavior of $\mathbf{w}$ and its feedback to the large-scale component, $\mathbf{v}$. With these new ingredients, the resulting ROM with the closure reads:

where $\dot{\mathbf{W}}_{v}$ and $\dot{\mathbf{W}}_{w}$ are the standard Gaussian white noises. The noise coefficients $\bm{\sigma}_{v}$ and $\bm{\sigma}_{w}$ are constant matrices and for simplicity they are diagonal. The four matrix or vector functions in the closure terms $\bm{\tau}^{(v)}_{w}(\mathbf{v})$, $\bm{\tau}^{(v)}_{v}(\mathbf{v})$, $\bm{\tau}^{(w)}_{w}(\mathbf{v})$, and $\bm{\tau}^{(w)}_{v}(\mathbf{v})$ can be arbitrarily nonlinear functions of $\mathbf{v}$. It is worthwhile to highlight that the conditional linear functions of $\mathbf{w}$ are not linear functions of $\mathbf{w}$. The conditional linearity means that $\mathbf{w}$ is linear with $\mathbf{v}$ being fixed. Therefore, the closure in (3) includes the effects from the nonlinear interactions between $\mathbf{v}$ and $\mathbf{w}$. Here, for the consistency with the quadratic nonlinearity in the starting model (2), the closure terms $\bm{\tau}^{(v)}_{w}(\mathbf{v})\mathbf{w}+\bm{\tau}^{(v)}_{v}(\mathbf{v})$ and $\bm{\tau}^{(w)}_{w}(\mathbf{v})\mathbf{w}+\bm{\tau}^{(w)}_{v}(\mathbf{v})$ are assumed to be quadratic, and therefore

where $\mathbf{D}^{(v)}_{v\tau}$ and $\mathbf{D}^{(w)}_{v\tau}$ are constant vectors, while $\mathbf{B}^{(\cdot)}_{\cdot}$, $\mathbf{C}^{(\cdot)}_{\cdot}$, and $\mathbf{D}^{(\cdot)}_{\cdot}$ are constant matrices.

It is worthwhile to highlight that many complex nonlinear systems have already fitted into the framework of (3). Some well-known classes of the models are [51, 54]:

• 1.

  Physics-constrained nonlinear stochastic models. Examples include the noisy versions of Lorenz models, low-order models of Charney-DeVore flows, and a paradigm model for topographic mean flow interaction.

• 2.

  Stochastically coupled reaction-diffusion models in neuroscience and ecology. Examples include stochastically coupled FitzHugh-Nagumo models and stochastically coupled SIR epidemic models.

• 3.

  Multi-scale models in turbulence, fluids and geophysical flows. Example include the Boussinesq equations with noise and stochastically forced rotating shallow water equation.

Other applications that share similar motivations or ideas as (3) include predicting the intermittent time-series of the monsoon intraseasonal variabilities [55, 56], recovering the turbulent ocean flows with noisy observations from Lagrangian tracers [57, 58], cheap solvable forecast models in dynamic stochastic superresolution [59, 60], and stochastic superparameterization for geophysical turbulence [8]. These examples further justify that the new ROM framework (3) is appropriate to characterize or approximate many nonlinear and non-Gaussian systems in various disciplines.

One of the major goals of building a ROM is to use it as a cheap forecast model. For turbulent systems, the forecast skill relies heavily on the accuracy of the initial values. However, it is often the case that only partial observations are available for complex systems. These partial observations usually correspond to the large-scale modes as they are the easiest to identify, while the variables that are below certain spatial scales are often hard to beapproximated. Therefore, it is important to combine the ROM with the partial observations for recovering the unobserved state variables in an accurate fashion,which then can be used as initializations for the forecast [61, 62, 63, 64, 65, 66, 67]. This is known as data assimilation (or filtering) [47, 68, 49, 46, 48].

Ensemble based data assimilation schemes [46, 69], such as the ensemble Kalman filter, have been widely utilized in practiceeliminate the rest of this sentence? and lead to many successes. However, there are a few intrinsic difficulties in applying these classical data assimilation methods. The main challenge comes from the sampling error as the state estimation is completely determined by the ensembles. To prevent the filter divergence, localization, noise inflation, and many other ad hoc tuning procedures are usually inevitable as effective remedies when using the ensemble based methods [70, 71]. In addition, if a large number of ensembles are used in the system to guarantee the accuracy of the data assimilation results, then the computational cost can significantly increase.

Different from the traditional data assimilation methods, the new ROM in (3) allows an efficient and accurate data assimilation scheme. The solution of the data assimilation result is given by closed analytic formulae and therefore it avoids all the ad hoc tuning procedures. In fact, the ROM in (3) can be written into a more abstract form,

where $\mathbf{A}_{0},\mathbf{a}_{0},\mathbf{A}_{1},\mathbf{a}_{1},\mathbf{B}_{1}$, and $\mathbf{b}_{2}$ are matrices and vectors that can depend nonlinearly on the state variables $\mathbf{v}$. Despite being highly nonlinear and non-Gaussian, the process of $\mathbf{w}$ in (5) is linear by design if a realization of the time series $\mathbf{v}(s)$ for $s\in[0,t]$ is given. Due to such a conditional linear property, the conditional distribution of $\mathbf{w}(t)$, given a realization of $\mathbf{v}(s)$ for $s\in[0,t]$, is Gaussian,

One desirable feature of the ROM framework (5) is that the conditional mean $\bm{\mu}$ and the conditional covariance $\mathbf{R}$ can be solved via the following closed analytic formulae [72, 51, 54]

The explicit formulae in (6)–(7) correspond to the optimal nonlinear filter (namely the online data assimilation) solution of the state variable $\mathbf{w}(t)$ given a realization of the observed time series $\mathbf{v}(s)$ for $s\in[0,t]$. Thus, $\bm{\mu}$ and $\mathbf{R}$ are also known as the filter posterior mean and the filter posterior covariance, respectively. The closed analytic formulae for the filtering solution allow an efficient and accurate online state estimation of the unobserved state $\mathbf{w}$ from the observed time series of $\mathbf{v}$ and the model structure of the ROM. The real-time recovery of the unobserved state provides the initial condition of $\mathbf{w}$, which together with the observed initial value of $\mathbf{v}$, advances the forecast utilizing the ROM. Note that the classical Kalman-Bucy filter [73] is the simplest special example of (7).

Due to the solvable conditional Gaussian statistics in (6)–(7), the new ROM is named as the conditional Gaussian ROM (CG-ROM).

What remains is to determine the matrices and vectors in the ROM (3), which correspond to the model parameters. As was discussed above, partial observation means only a time series of the observed large-scale variable $\mathbf{v}$ is available in practice. Yet, it is assumed that a time series of both $\mathbf{v}$ and $\mathbf{w}$ is accessible in the model calibration stage to simplify the procedure of determining the model parameters. This can be explained as using the reanalysis data in many practical applications, where the historical trajectories of the unobserved data has been inferred via certain Bayesian approaches [74]. Otherwise, an iteration process involving an alternation between state estimation and parameter estimation can be adopted for the model calibration. The basic algorithm of parameter estimation remains the same but the entire process is more complicated, which will diverge from the ROM development here, and is therefore treated as future work(see the discussion in Section 4). Nevertheless, only the value of the observed variable $\mathbf{v}$ is available in the forecast stage such that data assimilation becomes essential in recovering the initial value of $\mathbf{w}$.

The parameters in the ROM (3) can be divided into three sets:

1. (A).

   The parameters in the basic dynamics: $\mathbf{A}^{(v)}_{v}$, $\mathbf{B}^{(v)}_{vv}$, $\mathbf{B}^{(v)}_{vw}$, $\mathbf{A}^{(w)}_{w}$, $\mathbf{B}^{(w)}_{vv}$, $\mathbf{B}^{(w)}_{vw}$ in (3)

2. (B).

   The parameters in the deterministic part of the closure: $\mathbf{B}^{(\cdot)}_{\cdot}$, $\mathbf{C}^{(\cdot)}_{\cdot}$ and $\mathbf{D}^{(\cdot)}_{\cdot}$ in (4)

3. (C).

   The parameters in the stochastic part of the closure: $\bm{\sigma}_{v}$ and $\bm{\sigma}_{w}$ in (3)

The parameters in Set (A) are inherited from the FOM if the FOM is available. Otherwise, they can be handled in the same way as those in Sets (B) and (C), which are computed using the following optimization algorithm. For the simplicity of discussion, consider the discrete approximation of the original continuous version of the ROM (3) by applying the Euler-Maruyama scheme:

where $\varepsilon^{j}_{\mathbf{v}}$ and $\varepsilon^{j}_{\mathbf{w}}$ are standard independent and identically distributed Gaussian random variables at time $t_{j}$. Assume the noise coefficients are constants, namely $\bm{\sigma}_{v}^{j}:=\bm{\sigma}_{v}$ and $\bm{\sigma}_{w}^{j}:=\bm{\sigma}_{w}$ for all $j$. Denote $\bm{\Sigma}$ a block diagonal matrix with entries $(\bm{\sigma}_{v})^{2}\Delta t$ and $(\bm{\sigma}_{w})^{2}\Delta t$. Note that in many practical situations $\bm{\Sigma}$ is a diagonal matrix, from which $\bm{\sigma}_{v}$ and $\bm{\sigma}_{w}$ can be inferred. Further denote

where $\bm{\theta}_{L}$ contains all the parameters in $\bm{\tau}^{(v)}_{w},\bm{\tau}^{(v)}_{v}$ and $\bm{\theta}_{S}$ contains all those in $\bm{\tau}^{(w)}_{w},\bm{\tau}^{(w)}_{v}$. Here, the parameters are assumed to appear as multiplicative pre-factors of functions consisting of state variables $\mathbf{M}^{j}$ such that the closure terms can be written in an abstract form $\mathbf{M}^{j}\bm{\theta}$. The vector $\mathbf{C}^{j}$ in (9) includes all those terms on the right hand side of (3) that do not contain parameters. Then the optimization problem based on a linear regression can be formulated as [75]:

where $\bm{z}^{j+1}=[\mathbf{v}^{j+1};\mathbf{w}^{j+1}]$ and $J$ is the number of training snapshots. To solve the optimization problem (10), we use the iterative method to compute the following:

Physics constraints correspond to the conservation of the energy in the quadratic nonlinear terms in (1), namely, $\mathbf{u}\cdot B\left(\mathbf{u},\mathbf{u}\right)=0$. Incorporating the physics constraint into the ROM is central in preventing finite-time blow-up of the solutions and facilitates a skillful medium- to long-range forecast [38, 39, 50].

The physics constraints (or in general any constraint) for the model can be written as the following matrix equality:

where $\mathbf{H}$ contains those functions that are involved in the constraints. Then the optimization problem in (10) is reformulated as a constrained optimization problem,

The constrained optimization problem can be solved by the Lagrangian multiplier method [76], where the Lagrangian function yields

The solution can be obtained using the iterative method:

where

In Figure 1, we summarize all the models and corresponding equations used in this paper.

The first application of the CG-ROM framework is the following viscous stochastic Burgers equation [52, 53],

which is posed on the spatial domain $(0,L)$ supplemented with the Dirichlet boundary conditions $u(0,t)=u(L,t)=0$ for $t\geq 0$ and a suitable initial condition $u(x,0)=u_{0}(x)$ for $x\in(0,L)$. In (22), $\nu>0$ is the viscosity constant, $\lambda$ is the linear drag coefficient, $\gamma$ is the advection strength, $\sigma$ stands for the level of the stochastic noise, and $\dot{W}_{t}$ is a spatiotemporal white noise source.

The viscous stochastic Burgers equation (22) is very different from the standard deterministic inviscid Burgers equation, which is well-known for the development of discontinuities, namely shock waves. In (22), the external random forcing interacts with the deterministic nonlinearity to induce intermittent bursts of the signals in the form of viscous shocks, which will be dissipated via the viscous term before the appearance of discontinuities. Therefore, the viscous stochastic Burgers equation is a suitable test model that mimics turbulent behavior in nature with intermittency. It is also noticeable that the stochastic forcing is triggering turbulent features.

The second application is the quasi-geostrophic (QG) equation:

which is a widely used model that describes the large scale ocean circulation [1, 9]. In the QG equation (32), $\omega$ is the vorticity, $\psi$ is the streamfunction, $F$ is the external forcing, and theJacobian is defined as $J(A,B)=\partial A/\partial x\,\partial B/\partial y-\partial A/\partial y\,\partial B/\partial x$. There are two non-dimensional numbers in (32): $Re$ is the Reynolds number and $Ro$ is the Rossby number. The former is the ratio of inertial forces to viscous forces, representing the strength of turbulence, while the latter stands for the ratio of inertial force to Coriolis force. Note that, different from the viscous stochastic Burgers equation, the QG equation is a deterministic system and its turbulent features are induced completely by the nonlinear interactions between the state variables at different spatial scales.

In the following, a symmetric double-gyre wind forcing is imposed to the QG equation (32) [78, 79, 80, 81] with $F=\sin(\pi(y-1))$. A rectangular domain is considered here with $\Omega=[0,1]\times[0,2]$. Homogeneous Dirichlet boundary conditions are applied for both $\psi$ and $\omega$:

The two non-dimensional numbers are set to be $Re=450$ and $Ro=0.0036$, such that the system has moderate turbulent behavior and a fast rotation.