# UP-dROM - Uncertainty-Aware and Parametrised dynamic Reduced-Order Model – application to unsteady flows

Ismaël Zighed, Nicolas Thome, Patrick Gallinari, Taraneh Sayadi

► To cite this version:

Ismaël Zighed, Nicolas Thome, Patrick Gallinari, Taraneh Sayadi. UP-dROM - Uncertainty-Aware and Parametrised dynamic Reduced-Order Model – application to unsteady flows. Physical Review Fluids, 2025, 10 (11), pp.114902. 10.1103/f6ty-t6gl. hal-05416989

HAL Id: hal-05416989

https://hal.science/hal-05416989v1

Submitted on 16 Dec 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# UP-dROM - Uncertainty-Aware and Parametrised dynamic Reduced-Order Model – application to unsteady flows

Ismaël Zighed¹,², Nicolas Thome²,³, Patrick Gallinari², Taraneh Sayadi⁴

¹ Institut Jean le Rond d'Alembert - Sorbonne Université, ² ISIR - Sorbonne Université, ³ Institut Universitaire de France, ⁴ M2N - Conservatoire National des Arts et Métiers

Corresponding authors: Ismaël Zighed (ismael.zighed@sorbonne-universite.fr), Nicolas Thome (nicolas.thome@sorbonne-universite.fr), Patrick Gallinari (patrick.gallinari@sorbonne-universite.fr), Taraneh Sayadi (taraneh.sayadi@lecnam.net)

## Abstract

Reduced order models (ROMs) play a critical role in fluid mechanics by providing low-cost predictions, making them an attractive tool for engineering applications. However, for ROMs to be widely applicable, they must not only generalise well across different regimes, but also provide a measure of confidence in their predictions. While recent data-driven approaches have begun to address nonlinear reduction techniques to improve predictions in transient environments, challenges remain in terms of robustness and parametrisation.

In this work, we present a nonlinear reduction strategy specifically designed for transient flows that incorporates parametrisation and uncertainty quantification. Our reduction strategy features a variational autoencoder (VAE) that uses variational inference for confidence measurement. We use a latent space transformer that incorporates recent advances in attention mechanisms to predict dynamical systems. Attention's versatility in learning sequences and capturing their dependence on external parameters enhances generalisation across a wide range of dynamics. Prediction, coupled with confidence, enables more informed decision-making and addresses the need for more robust models.

In addition, this confidence is used to cost-effectively sample the parameter space, improving model performance a priori across the entire parameter space without requiring evaluation data for the entire domain.

**Keywords:** Reduced Order Modelling, Dynamical Systems, Deep Learning, Uncertainty Quantification

## 1. Introduction

The integration of machine learning with physics has fostered a mutually beneficial relationship, driving significant advancements in computational modelling. This synergy has given rise to the development of both physics-informed and physics-inspired models. Physics-informed models, like PINNs [1], explicitly incorporate physical laws into their structure, ensuring that solutions adhere to known scientific principles [2]. On the other hand, physics-inspired models, such as Hamiltonian and Lagrangian networks, are based on the conceptual framework of classical mechanics and other physical theories. These models draw inspiration from the rich mathematical structures found in physics, using concepts like energy conservation, action minimisation, and symmetries to inform their architecture and training procedures. By leveraging these fundamental principles, physics-inspired models are able to encode essential physical behaviours into their design, often leading to enhanced generalisation and more physically plausible predictions. For example, Hamiltonian networks use the principles of Hamiltonian mechanics to model the dynamics of a system, while
Lagrangian networks leverage the idea of least action to optimise their predictions [3, 4]. These approaches allow for the modelling of complex physical systems by inherently incorporating dynamical properties and constraints, without requiring explicit equations for every interaction, thus capturing the essence of physical systems in a computationally efficient manner.

Building on the integration of machine learning and physics, the concept of Neural Ordinary Differential Equations (Neural ODEs) has recently emerged as a powerful tool for modelling dynamic systems [5]. Neural ODEs provide a continuous-time framework for representing the evolution of dynamic systems, where the time evolution of states is modelled by neural networks in a way that is both smooth and adaptive. This continuous formulation contrasts sharply with traditional discrete-time architectures, such as Residual Networks (ResNets) [6], which operate on fixed time steps. Extending these ideas, Neural Operators aim to learn mappings between function spaces rather than finite-dimensional vectors [7, 8, 9]. This innovative approach enables the learning of operators that map between functions, making it highly effective for solving partial differential equations (PDEs) that describe complex systems across different domains.

One notable application of these advanced modelling techniques is Reduced Order Modelling (ROM). ROMs aim to offer computationally efficient alternatives to large, data-intensive models, which is especially valuable for real-time simulations and control in critical systems [10]. By simplifying the system representation, ROMs improve interpretability and transparency, key factors in applications where understanding the underlying mechanics is just as crucial as making accurate predictions. Traditional ROM approaches, such as Proper Orthogonal Decomposition (POD) [11], achieve dimensionality reduction by projecting data onto a lower-dimensional manifold based on the overall contribution to the energy norm. An alternative to POD, Dynamic Mode Decomposition (DMD) [12], focuses on decomposing complex systems into the most dominant time-evolving modes. While both POD and DMD, along with their variations [13, 14], have proven effective in identifying coherent structures within data, they are limited by their linear nature, which restricts their ability to fully capture nonlinear dynamics. Consequently, there has been a growing shift towards nonlinear methods [15, 16]. In this context, machine learning techniques have been increasingly used to both identify the reduced basis [17] and to represent the dynamics within the lower-dimensional manifold [18]. A recent advancement in this area is the machine-learning closure model for POD, referred to as CD-ROM [19], where the projection error orthogonal to the POD basis is addressed by a historical term predicted through a Neural ODE. This methodology offers a more accurate representation of the system's dynamics, significantly enhancing the performance of reduced order models in capturing the complexities of nonlinear systems.

Koopman theory provides a robust mathematical framework for dealing with complex dynamical systems by assuming the existence of an infinite-dimensional space in which the dynamics are linear. This concept has been applied in various machine learning implementations [20], but it requires increasing the dimensionality towards infinity, which is often impractical and different from the Reduced Order Model (ROM) framework that our work seeks to integrate. Work combining ROMs with Koopman theory has shown impressive results on complex dynamical problems such as the Navier-Stokes equations [21, 22].

However, these methods have two major limitations. First, they need to be parametrised, or in other words learn and generalise across different operating conditions, referred to in this paper as 'parameters'. Second, they are inherently deterministic, which limits their ability to account for uncertainty in predictions. To address these limitations, we emphasise the importance of adopting a probabilistic framework. Such a framework allows for the quantification of uncertainty, which is critical for providing confidence levels in predictions. A probabilistic approach also offers improved robustness to model inaccuracies and support for risk-aware decision-making. These benefits are critical for fostering confidence and enabling the integration of machine learning models into model-based design.

Recent advancements have focused on integrating probabilistic methods with Reduced Order Models (ROMs) using techniques like Laplace Approximation [23], Variational Inference [24, 25], Gaussian Processes [26], Ensemble Methods [27], and Monte-Carlo Dropout [28]. However, these methods are seldom applied to non-stationary dynamical systems in practice, largely due to the challenges of integrating them with temporal models such as Neural Ordinary Differential Equations (NODEs), Recurrent Neural Networks (RNNs), or Transformers. An exception is Gaussian Processes, which can model time-series data as multivariate Gaussian distributions, but they come with high computational costs and often require complex tuning of
kernel parameters.

To address the difficulty of capturing time dependencies, several approaches have combined traditional reduction methods with Transformers in the latent space [29]. Transformers, in particular, have proven highly effective for parametrisation via cross-attention [30], allowing models to adapt dynamically to different excitation regimes driven by external variables. Additionally, Transformers have shown their versatility and robustness beyond the large language model (LLM) framework, emerging as the state-of-the-art in time-series forecasting. Their generality and flexibility have also facilitated their growing adoption in partial differential equation (PDE) applications, demonstrating their potential to manage and model complex dynamical systems.

Through UP-dROM, we aim to contribute to these advances by proposing a Transformer-based, dynamic, parametric and variational ROM. This model is designed to solve unsteady PDEs while considering and adapting to external excitation variables through cross-attention mechanisms, thus exploiting breakthroughs in large language models (LLMs) for PDE applications. In addition, it provides a probabilistic counterpart to most deterministic dynamic modelling baselines, enabling uncertainty quantification. Crucially, our approach achieves this in a computationally efficient manner, significantly reducing the cost compared to Gaussian processes or Bayesian models, while maintaining robustness and scalability for complex dynamical systems.

The paper is structured as follows: Section 2 presents the various steps involved in constructing the model, beginning with the reduction procedure and the corresponding dynamical representation within the reduced manifold. This section also includes the model validation using the test case of a laminar, unsteady flow around a bluff body. The process of parametrising the model is then outlined in Section 3. Section 4 describes the procedure for incorporating uncertainty measures into the model, considering both applications in physical space and parameter space. Section 5 describes the adaptive procedure, which uses the uncertainty measure to efficiently sample the parameter space. Finally, Section 6 concludes the paper with a discussion on the findings and directions for future work.

## 2. Components of UP-dROM

Consider a dynamical system represented by the equation:

$$ \frac{d\phi}{dt} = F(\phi, \xi, t), \qquad (1) $$

where $F$ is a nonlinear operator governing the spatial and temporal evolution of the system, $\phi$ represents the state variables and $\xi$ the parameters governing the dynamical behaviour of the system. The objective is to develop a reduced-order model that accurately represents the dynamics of the given partial differential equation. Considering the high dimensionality of the problem and the inherent presence of coherent structures (modes) in most dynamical systems, we employ dimensionality reduction techniques to both preserve the coherence and reduce the system's dimensionality. To effectively analyse and capture these space-time coherent structures, we propose a space-time strategy that compresses spatial information while simultaneously modelling the temporal dependencies within this lower-dimensional manifold.

### 2.1. Dimensionality reduction

In this work, we reduce the spatial component of the system using a Variational Autoencoder (VAE), which facilitates the compression of high-dimensional physical data into a lower-dimensional nonlinear manifold. The choice of a VAE over a traditional autoencoder is motivated by its ability to incorporate uncertainty quantification, which is a critical aspect of our approach, as discussed in Section 4. The VAE is particularly advantageous because it introduces a probabilistic framework, enabling the model to not only capture the underlying structure of the data but also quantify the uncertainty in the learned representations [24]. This probabilistic nature allows for more reliable inference, especially in situations where uncertainty plays a significant role in prediction or system behaviour. Variational Autoencoders offer the advantage of learning continuous, disentangled latent representations. This structured latent space improves the robustness of
inference and enhances generalisation beyond the training distribution, especially when compared to traditional Auto-Encoders. The ability to generalise more effectively is a key factor driving the widespread use of VAEs in generative AI frameworks [31].

Moreover, VAEs and Autoencoders (AEs) have been shown to outperform traditional methods in compressing data from complex physical systems [21, 24, 32, 33, 25]. This is because VAEs and AEs do not impose strict orthogonality constraints, which is a key limitation in traditional linear compression techniques, such as Proper Orthogonal Decomposition (POD), PCA [11, 34]. In contrast, VAEs and AEs offer the flexibility to relax this constraint, making them better suited for capturing intricate nonlinear structures within the data. Furthermore, the absence of a linearity constraint enables VAEs to generalise more effectively to unseen data distributions. This ability to generalise is particularly important when dealing with real-world physical systems, where the underlying dynamics may vary across different regimes or boundary conditions.

Overall, the use of VAEs provides significant advantages in terms of both data compression and generalisation performance, especially in nonlinear, high-dimensional settings. Their probabilistic nature allows uncertainty quantification, making them a more powerful tool compared to traditional methods that rely on linear approximations. This robustness and flexibility make VAEs a promising approach for modelling and analysing complex physical systems, where the capturing of nonlinear dynamics and the ability to handle uncertainty are crucial for accurate predictions and decision-making. The architecture of the model is shown in Figure 1. In the provided architecture, $\mu$ and $\sigma$ represent the mean and standard deviation of the latent variable distribution, respectively. These parameters define the probabilistic mapping from the input space to the latent space. The determination of each parameter is further explained in Section 4.

## 2.2. Temporal evolution

To capture the time component of these structures within the latent space, several strategies can be employed. These include utilizing a Koopman operator [22], solving the governing equations analytically using Galerkin projection in conjunction with a deep-learning closure model [19], or applying Long Short-Term Memory (LSTM) networks to capture both short-term and long-term temporal dependencies [35]. In this study, we propose the use of Transformers, similar to those employed in large language models (LLMs) [30]. This choice is driven by two primary considerations: (i) the attention mechanism in Transformers serves as an effective memory component, facilitating the model's ability to retain relevant information over time, and (ii) the combination of self-attention and cross-attention mechanisms enables the model to capture complex and nonlinear dependencies, including interactions with external variables. This approach enhances the model's generalisation capabilities and overall robustness. The resulting architecture is illustrated in Figure 1.

Regarding point (i), The limitations of simple Markovian models have been clearly demonstrated by the Mori-Zwanzig formalism [36]. Specifically, in a test case similar to that considered in this paper, Gupta et al. [22] showed that purely Markovian models fail to capture the finite-dimensional invariant Koopman subspace necessary to fully describe the system dynamics. Their results emphasise the importance of incorporating the effects of unresolved observables on the finite resolved observable/latent space via a memory term. In other words, when constructing a low-dimensional latent space, memory terms are necessary to capture the full dynamics of the dynamical systems. This is also related to Takens' theorem, [37] which proves that the state of a dynamical system can be reconstructed from a set of time-delayed embeddings based on limited observations. The projection performed by the VAE maps the high-dimensional observations $\Phi$ onto a lower-dimensional latent space $z$. Time-delayed embeddings constructed from the latent representation $z$ thus serve as sparse observations of the full state $\phi$, allowing the application of Takens' theorem directly within the latent space. The trajectory of the system can be uniquely reconstructed from the delayed observations:

$$z(t), z(t - \tau), z(t - 2\tau), \dots, z(t - (d - 1)\tau) \qquad (2)$$

By using Takens' theorem, we ensure that the sequence of latent vectors encapsulates the full state of the dynamical system. The transformer, with its attention mechanism, models the interactions and dependencies within this sequence, allowing accurate prediction of future system states and capturing complex temporal
<figure><img src="image_2.png" /><figcaption>Figure 1: UP-dROM architecture with q the look-back window</figcaption></figure>

dynamics. Similar to how Large Language Models (LLMs) use attention to capture long-range dependencies and contextual relationships in sequential data, transformers in this context exploit time-dependent structures in latent space. Transformers are often considered superior for handling time sequences and sequence self-dependencies compared to state-of-the-art Markovian or recurrent networks such as LSTMs. Unlike recurrent architectures, which process data sequentially and thus suffer from vanishing gradient issues or limited memory, Transformers use attention mechanisms to model dependencies regardless of their position in the input sequence [30]. Transformers can capture both the Markovian dynamics and the memory effects effects, as well as external dependencies through cross-attention, within a unified framework, making them well-suited for reconstructing the dynamics in the latent space. The VAE and the transformer are trained jointly, resulting in the model's loss function, shown in equation 3. The full pass loss refers to direct supervision using a state vector that is encoded, propagated through time, and decoded. The VAE's loss function is described in Section 4, while the overall model's loss is detailed in Appendix 7, along with the relevant hyperparameters. Our complete implementation is available on Github [38].

$$ \mathcal{L} = \mathcal{L}_{VAE} + \mathcal{L}_{transformer} + \mathcal{L}_{full\ pass} \qquad (3) $$

## 2.3. Validation of UP-dROM – application to flow around a cylinder

The model is validated by predicting the evolution of flow around a blunt object, specifically a cylinder. The flow dynamics are governed by the incompressible Navier-Stokes equations:

$$ \nabla \cdot u = 0 \qquad (4) $$

$$ \frac{\partial u}{\partial t} = -\nabla p - (u \cdot \nabla)u + \frac{1}{Re} \nabla^2 u \qquad (5) $$

Where u is the velocity field, p is the pressure, and Re is the Reynolds number. The data are generated using a solver based on the Immersed Boundary Method [39]. This solver utilises two grids: a regular Eulerian grid that spans the entire flow domain and a Lagrangian grid that conforms to the boundary of the bluff body, enabling accurate representation of the obstacle-fluid interaction. Once generated, only the Eulerian grid data of the flow field are used; the model is not explicitly informed of the bluff body's position or shape and must learn it implicitly from the data.

The dataset provided to the model publicly available in [40], includes all Eulerian grid points and starts from the stationary unstable solution and captures the transient dynamics leading to the stable limit cycle
every time increment $dt$. It includes velocities $u$ and $v$ in the $x$ and $y$ directions, respectively, stacked in a state vector $\Phi$. The dataset is divided into a training and a test sets. In this example, we split the data in half by assigning even time steps to the training set and odd time steps to the test set, ensuring equal representation of the transients and the limit cycle in both the training and test sets. This results in an equal split of 50% between the training and test data, while preserving the number of oscillation periods and balancing the presence of transients and limit cycles in each set. Consequently, the effective time interval between consecutive data points in either set doubles to $2dt$. The model is trained using fixed-size look-back windows and learns to predict multiple future steps in an autoregressive manner as illustrated by Figure 1. During inference, starting only from the initial condition and the unstable steady state, the model generates the entire dynamical trajectory by rolling out predictions step by step. In this test case, the VAE projects the velocity fields onto a four-dimensional latent space and the Transformer predicts the trajectory, as shown in figure 2. This figure illustrates the evolution from the initial condition (IC), which corresponds to an unstable fixed point acting as a repeller. This point gives rise to a low dimensional manifold bounded by the characteristic periodic orbit of the cylinder wake flow, as represented in the identified latent space.

<figure><img src="image_3.png" /><figcaption>Figure 2: Latent manifold visualization in three different 3D projections with the 4th dimension as a colour gradient.</figcaption></figure>

The flow behaviour is primarily influenced by the Reynolds number $R_e$, which is treated here as the external variable or parameter $\xi$. The Reynolds number is a dimensionless quantity that characterises the ratio of inertial forces to viscous forces in a fluid flow and is defined as:

$$R_e = \frac{\rho UL}{\mu}, \qquad (6)$$

where $\rho$ is the fluid density, $U$ is the characteristic velocity, $L$ is the characteristic length (such as the length of the bluff body), and $\mu$ is the dynamic viscosity of the fluid. For this particular bluff body, when $R_e$ remains below its critical value, the flow closely resembles a stationary solution, with minimal deviations over time. However, beyond this critical Reynolds number, vortex shedding and unsteady flow dynamics emerge in the wake of the cylinder, leading to a more complex and time-dependent behaviour. Figure 3 shows the base flow provided to the model at the initial time $t_0$ and a representative frame at a later time $t$ when the limit cycle is established. The model should be able to capture the transition dynamics and the stable limit cycle dynamics.
<figure><img src="image_4.png" /><figcaption>Figure 3: Predicted kinetic energy signal using UP-dROM compared with the ground truth.</figcaption></figure>

The transition from the base flow to the stable limit cycle is captured through the evolution of the kinetic energy signal. In Figure 3, we compare the inferred kinetic energy to the ground truth, followed by a forecasting window that demonstrates the model’s performance. The kinetic energy can be computed from the state vector $\phi$, which includes both the velocity fields $u$ and $v$ over time $t$, with the total number of points in each direction given by $n_{xy}$. The formula for kinetic energy is given by:

$$k_t = \frac{1}{2n_{xy}} \sum_{d=1}^{n_{xy}} (u_{d,t}^2 + v_{d,t}^2) \qquad (7)$$

This metric effectively averages the spatial components in a way that is physically meaningful. By eliminating the spatial dimensions, it enables us to concentrate on time-dependent signals, making it a valuable visualisation tool for this study. The results show that the model performs comparably to similar architectures and strategies that use delay embedding to represent the system dynamics in the latent space [22, 19].

A key decision in our approach is the preference for a Transformer model over other time-embedding architectures, such as Long Short-Term Memory (LSTM) networks. This choice is informed by recent advancements in embedding memory mechanisms, which significantly enhance the model’s predictive capabilities [29, 41, 25, 42]. Specifically, Transformers leverage the attention mechanism, which allows the model to capture long-range dependencies in the time-series data more effectively than traditional methods. In addition to self-attention, the incorporation of cross-attention mechanisms in Transformers enables the model to respond adaptively to external variables or influences, providing a more comprehensive understanding of the system dynamics. This architecture choice also supports the broader objective of improving model generalisation. By enabling the model to learn relationships and interactions between the data and an external variable such as $R_e$, Transformers are particularly well-suited to address the challenges of generalizing across different conditions, thereby enhancing the robustness and versatility of the model in diverse physical system scenarios.

### 3. parametrisation

One of the primary applications of reduced-order models (ROMs) is in multi-query problems, such as control, where a function is evaluated across a wide range of operating conditions to identify the optimal operating point in accordance with a specific design objective. ROMs are particularly advantageous in this context due to their ability to provide cost-effective function evaluations, making them the preferred choice for conducting such analyses. However, for ROMs to be truly effective in these applications, they
must maintain predictive accuracy across a range of varying parameters. This requirement represents a significant challenge, as many data-driven reduced-order models applied to dynamical systems fail to satisfy this criterion.

Consider the example of flow around a bluff body, as shown in Figure 3, which serves as a simple dynamical system with a varying Reynolds number, as depicted in Figure 4. As the Reynolds number changes, the system undergoes a bifurcation and transitions into an unstable state. This is evidenced by the steady-state branch bifurcating into oscillatory branches, with unsteady von Kármán vortex street patterns emerging in the flow domain. For the case of 2D flow around a bluff body considered here, we considered an elliptical obstacle with an aspect-ratio $\beta = 1.35$. This aspect-ratio sets the bifurcation at a critical Reynolds Number $R_{ec} \approx 60$.

<figure><img src="image_5.png" /><figcaption>Figure 4: Schematic of a bifurcation diagram with respect to the varying Reynolds number. The two inserted plots are the kinetic energy signals of two simulations performed at $Re = 50$ and $Re = 90$, before and after the critical Reynolds number, respectively.</figcaption></figure>

The introduction of a parameter therefore modifies the unstable behaviour of the system, as indicated in the Figure, by eliminating the vortex street, which in turn alters key flow properties such as the drag coefficient. For a reduced-order model to be effective in predicting such a dynamical system, it must be generalizable across different parametric regimes, which often exhibit significantly distinct dynamical behaviours.

The primary objective of this study is to parametrise the model introduced in the previous section and subsequently assess its predictive capabilities across a range of parameter values. To achieve this, we first introduce a parametrisation, denoted as $\xi$, to provide the model with information about external variables that may induce different dynamic response regimes. We incorporate the cross-attention mechanism within the Transformer architecture to account for these external variables during both the learning and inference processes. Specifically, each attention block includes a self-attention component to capture intra-sequence correlations, as well as a cross-attention component to integrate the external parameters. Moreover, to identify an appropriate latent representation for the various dynamic regimes, both the encoder and decoder are informed of these external variables, ensuring robust encoding and decoding across the entire parameter space similar to the approach proposed in [24]. The modification to the model architecture is illustrated in Figure 5.
<figure><img src="image_6.png" /><figcaption>Figure 5: UP-dROM parametrised by external excitation variable $\xi$.</figcaption></figure>

To evaluate the model’s performance, it was first trained using two Reynolds numbers, 90 and 120, both of which correspond to the post-transition limit cycle regime. The model’s performance at these two training points is highlighted in Figure 6, where the evolution of the total kinetic energy is compared. The Figure demonstrates that the introduction of the parameter does not affect the predictions of the parametrised model, which remain consistent with the performance observed in the previous section. This illustrates the effectiveness of the parametrisation in maintaining accurate inferences.

Additionally, Figure 6 shows the model’s performance when applied to a regime prior to the critical Reynolds number, a regime on which it was not trained. In this regime, the dynamics are significantly different, corresponding to the stable fixed point where the flow is steady and no limit cycle behaviour is observed. This is evident from the evolution of the kinetic energy, which exhibits a decaying and stable trend.

Remarkably, the model captures the underlying dynamics even though it was trained on only a subset of the possible behaviours. For example, it correctly identifies the transition region where the behaviour shifts to a stationary-like state. At first glance, it is surprising that the model can accurately capture bifurcation behaviour and infer the critical Reynolds number. This marks the point at which the real part of the eigenvalues of the Jacobian, evaluated at the fixed point, crosses zero. This indicates the onset of instability and the transition from the fixed point to a periodic orbit. Previous studies, such as Strykowski & Sreenivasan [43], demonstrated a linear relationship between the Reynolds number and the instability growth rate, which was later confirmed numerically by Canuto & Taira [44]. As this relationship is linear, only two data points are needed to estimate the growth rate mapping with the Reynolds number. We believe that this simple relationship, coupled with the fact that the model is parametrised and learns across the parametric phase-space, allows the model to extrapolate from just two training points and correctly infer the dynamical state of the system for out-of-distribution Reynolds numbers, such as Re = 50. Therefore, while the overall energy level is not accurately predicted, the model successfully captures the correct dynamical behaviour. In the following section, we propose a strategy to adaptively correct the minor errors in the model’s predictions.
<figure><img src="image_7.png" /><figcaption>Figure 6: Inferred and true kinetic energy signals in the in-distribution and out-of-distribution training data, superimposed on a sketch of the bifurcation plot with varying Reynolds number.</figcaption></figure>

To test whether the performance of UP-dROM extends to a larger dimensional parameter space, the aspect ratio of the bluff body was added as a second parameter. The aspect-ratio, as defined in Figure 7, characterises the shape of the bluff body and significantly influences the flow behaviour. An aspect-ratio smaller than 1 indicates that the obstacle is taller than wide, creating a stronger disruption in the flow. This increased blockage leads to earlier vortex shedding and a bifurcation occurring at a lower Reynolds number. In contrast, an aspect-ratio greater than 1 corresponds to an elongated body, similar to an airfoil, which disturbs the flow less significantly. In this case, the onset of vortex shedding and bifurcation is delayed, occurring at a higher Reynolds number due to the more streamlined geometry. This effect is illustrated in Figure 7, where different aspect ratios result in distinct flow regimes, even though the Reynolds number remains unchanged.

<figure><img src="image_8.png" /><figcaption>Figure 7: Aspect ratio definition (a) and $u$ velocity field at $R_e = 50$ for aspect-ratio $\beta = 1.5$ (b), and $\beta = 0.75$ (c).</figcaption></figure>

Performance with a two-variable parametrisation ($R_e, \beta$) can be visualised using a pointwise scaled MSE, as shown in Figure 8, which ensures that each grid point is equally weighted in the evaluation. The scaled MSE, computed on test data, is obtained locally using the magnitude of local variations. The resulting
heatmap provides a *a-posteriori* perspective on the model’s ability to generalise and its overall performance. The signals shown on the heatmap represent the evolution of the kinetic energy for that configuration over time. The model performs well on the training sets, demonstrating its ability to learn and predict. In particular, its interpolation results are strong. However, similar behaviour to the single parameter case is observed. UP-dROM is able to distinguish between stable and limit-cycle dynamics in the two-parameter regime. However, the accuracy of the prediction deteriorates when the model is used in extrapolation mode.

<figure><img src="image_9.png" /><figcaption>Figure 8: Scaled MSE in the parameter space with predicted (blue) vs true (red) kinetic energy signals at various location of the ($R_e, \beta$) parameter space.</figcaption></figure>

# 4. Uncertainty quantification

Most conventional reduction strategies rely on deterministic estimation, where the latent space projection is represented by a single point estimate. This approach fails to capture the inherent uncertainty within the model. In contrast, UP-dROM utilises variational methods to assess uncertainties, providing a probabilistic projection in the latent space instead of a deterministic one. This is achieved through the use of a Variational Autoencoder (VAE). Specifically, given the observed data $X = D$ (i.e., the dataset $D$), the method seeks to infer the posterior distribution $p(z|X = D)$, where $z$ represents the latent variables.

The central idea connecting variational inference and uncertainty quantification revolves around approximating the true posterior distribution $p(z|X = D)$ using a surrogate distribution $q(z|X = D)$. In this case, the surrogate is chosen to be a Gaussian distribution due to its desirable mathematical properties and ease of interpretation. To ensure that the surrogate closely approximates the true posterior, we minimise the
Kullback-Leibler Divergence (KLD) between $p(z|X = D)$ and $q(z|X = D) \sim \mathcal{N}(\mu, \sigma^2)$, as:

$$ \text{KLD}(p(z|X = D) \parallel q(z|X = D)) = \int p(z|X = D) \log \frac{p(z|X = D)}{q(z|X = D)} dz. $$

Since $p(z|X = D)$ is not known a priori, we instead utilise the Evidence Lower Bound (ELBO), which can be computed at any stage. Variational inference theory states that minimizing the ELBO, also minimises the KLD between the true posterior and the surrogate distribution [31]. The ELBO is given by:

$$ \text{ELBO} = \text{KLD}(p(z) \parallel q(z|X)) + \mathbb{E}_{z \sim q}[\log(p(X|z))], \quad (8) $$

where $p(z)$ is the prior Gaussian distribution representing the distribution of the latent space. It is typically chosen as $p(z) \sim \mathcal{N}(0, I)$, which imposes a standard normal prior on the latent variables. This choice ensures that the latent space is continuous and well-structured. It also simplifies the mathematical treatment of the model, making the optimisation more tractable. The KLD between two Gaussian distributions such as $q(z|X = D)$ and $p(z)$ can be computed as follows,

$$ \text{KLD}(q(z|X = D) \parallel p(z)) = -\frac{1}{2} \sum_{i=1}^{d} (\sigma_i^2 + \mu_i^2 - 1 - \log(\sigma_i^2)), \quad (9) $$

$$ \text{if} \quad \begin{cases} q(z|X = D) \sim \mathcal{N}(\mu, \sigma^2) \\ p(z) \sim \mathcal{N}(0, I) \end{cases} $$

The term, $\mathbb{E}_{z \sim q}[\log(p(X|z))]$, corresponds to the expected log-likelihood of the data, which measures how well the decoder is able to reconstruct the input $X$ from the latent variable $z$. Assume the decoder likelihood is a multivariate Gaussian with diagonal covariance:

$$ p(\mathbf{x}|\mathbf{z}) = \mathcal{N}(\mathbf{X}; \hat{\mathbf{X}}, \sigma^2\mathbf{I}) $$

where $\hat{\mathbf{X}} = \mathcal{D}(z)$ is the decoder output and $\sigma^2\mathbf{I}$ is the covariance matrix.

$$ p(\mathbf{X}|\mathbf{z}) = \frac{1}{(2\pi)^{D/2}|\sigma^2\mathbf{I}|^{1/2}} \exp \left( -\frac{1}{2}(\mathbf{X} - \hat{\mathbf{X}})^{\top}(\sigma^2\mathbf{I})^{-1}(\mathbf{X} - \hat{\mathbf{X}}) \right) $$

$$ p(\mathbf{X}|\mathbf{z}) = \frac{1}{(2\pi\sigma^2)^{D/2}} \exp \left( -\frac{1}{2\sigma^2} \|\mathbf{X} - \hat{\mathbf{X}}\|^2 \right) $$

$$ \log p(\mathbf{X}|\mathbf{z}) = -\frac{D}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\|\mathbf{X} - \hat{\mathbf{X}}\|^2 $$

$$ \mathbb{E}_{z \sim q}[\log p(\mathbf{X}|\mathbf{z})] = -\frac{D}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\mathbb{E}_{z \sim q}[\|\mathbf{X} - \hat{\mathbf{X}}\|^2] $$

Maximizing this expected log likelihood is equivalent to minimizing the mean squared error (MSE):

$$ \mathbb{E}_{z \sim q}[\log(p(X|z))] = \text{MSE}(\hat{\mathbf{X}}, X), $$

A Variational Autoencoder (VAE) integrates the Evidence Lower Bound (ELBO) into its loss function, as

$$ \mathcal{L}_{VAE} = \text{ELBO}. \quad (10) $$

Therefore, the loss function of the VAE can be computed as

$$ \mathcal{L}_{VAE} = \text{MSE}(\hat{\mathbf{X}}, X) - \frac{1}{2} \sum_{i=1}^{d} (\sigma_i^2 + \mu_i^2 - 1 - \log(\sigma_i^2)). \quad (11) $$
Through this loss function, the model learns how to better approximate the true posterior distribution using the Gaussian surrogate parametrised by a mean $\mu$ and a variance $\sigma^2$. At the end of the optimisation process, $q(z|X = D)$ closely approximates the true posterior and provides the density function of the latent projection. Specifically, if the model is confident, $q(z|X = D)$ will be a narrow Gaussian distribution centred around the expected latent coordinates, indicating low uncertainty in the latent representation. Conversely, if the model is uncertain, $q(z|X = D)$ will have a wide standard deviation, reflecting higher entropy in the latent space. This broader distribution captures the range of possible latent coordinates consistent with the observed data, effectively quantifying the uncertainty in the projection.

During training, the encoder is used at each time step to continuously learn the best-fitting projection of the dynamics at each step. This is not the case during inference, where the dynamics are propagated solely within the latent space using auto-regression. After the initial encoding, the encoder and its associated $\sigma$ tensor, which captures the uncertainty, do not need to be recomputed, keeping the computational cost low throughout the prediction.

<figure><img src="image_10.png" /><figcaption>Figure 9: UP-dROM Uncertainty Quantification process after inference</figcaption></figure>

We leverage the probabilistic latent projection to compute the uncertainty of the predictions. To do so, a second pass in the model is introduced. In this pass, the predicted dynamic window is processed through the VAE as shown in Figure 9. At each time step, the encoder generates a latent distribution. The method then quantifies the uncertainty by measuring the dispersion $\nu$ of the output sampled $n$ times from this latent distribution. These samples are indexed by $i$ in Figure 9. This strategy provides a confidence measure for the entire prediction window, not just the initial state. It also separates the uncertainty estimation from the main forecasting process, allowing for a more comprehensive analysis without increasing the computational burden of the forecasting task. As a result, the strategy remains both robust and scalable, even for large-scale applications.

Uncertainty quantification (UQ) using this approach can therefore be interpreted as a stochastic process, where uncertainty is measured by the standard deviation $\nu$ of the output trajectories over multiple samples from the latent distribution. Mathematically, the uncertainty $\nu$ at a spatial location $d$, time $t$, and parametrisation $\xi$, given $n$ samples, is defined as:

$$ \nu_{d,t,\xi} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (\phi_{i,d,t,\xi} - \bar{\phi}_{d,t,\xi})^2} $$

where $\bar{\phi}_{d,t,\xi}$ is the ensemble mean. The measure of uncertainty depends on the spatial dimension $d$, the time
$t$ and the parametrisation $\xi$.

### 4.1. Uncertainty in space

We can use the uncertainty measure $\nu_{d,t,\xi}$ to visualise how uncertainty is distributed across the computational domain at different time steps for a fixed parametrisation. Figure 10-a shows the $u$-velocity fields at three key time steps: (1) before the wake forms, during the unstable stationary-like behaviour, (2) during the transition, and (3) once the system has reached the stable limit-cycle, as highlighted in the kinetic energy plot. Figure 10-b illustrates the distribution of uncertainty in space $\nu_d$ at these same time steps, computed as:

$$ \nu_{d,t=t_j} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (\phi_{i,d,t=t_j} - \bar{\phi}_{d,t=t_j})^2} $$

<figure><img src="image_11.png" /><figcaption>Figure 10: $U$ velocity field at the selected snapshots (a), the corresponding uncertainty fields $\nu_d$ (b), and the snapshot locations indicated on the kinetic energy signal (c), for a fixed parametrisation $\xi$: $R_e = 90$.</figcaption></figure>

This figure illustrates a strong correlation between the known space-time dynamics and the uncertainty measurement. Prior to the system's transitioning to the stable limit cycle, uncertainty is evenly distributed along the symmetry axis of the flow domain, as expected, but this distribution breaks down as the dynamics transition. It is in this transient regime that the uncertainty across the entire domain is the highest. The highest confidence is observed upstream of the obstacle, where the dynamics remain relatively stationary throughout the time window. In contrast, lower confidence is concentrated at the edges of the obstacle's wake, suggesting that uncertainty is more closely related to fluctuations in the flow, rather than the instantaneous flow values in specific regions. These fluctuations increase variability within the system, leading to higher uncertainty.

### 4.2. Uncertainty in Time

The uncertainty measurement $\nu_{t,d,\xi}$ can be used to design confidence intervals. By fixing the location at a particular point and a specific parametrisation, Figure 11 illustrates such confidence intervals on both velocity fields.
<figure><img src="image_12.png" /><figcaption>Figure 11: Confidence intervals for velocity fields, at a particular point *d* and a specific parametrisation ξ : R_e = 90.</figcaption></figure>

Notably, the true signal lies almost entirely within the confidence interval. It is important to emphasise that this Uncertainty Quantification (UQ) strategy is not time-dependent and is solely driven by the input state. If it were time-dependent, the uncertainty would gradually increase as the dynamic process evolves autoregressively. Instead, we propose a likelihood measure for obtaining a given value, independent of previous values. This approach provides an instantaneous uncertainty measure with minimal computational cost.

As shown in Figure 11, the greatest uncertainties occur at the peaks of the signal, corresponding to instances when the local dynamics are most critical. This typically happens when the obstacle's wake shifts within the flow domain, causing substantial changes that the model must predict with accuracy. These transitions trigger increased uncertainty. This behaviour has a significant physical interpretation: uncertainty is highest in regions and time periods where the flow dynamics are most sensitive to changes, such as during sharp transitions, shifts in wake structure, small-scale patterns, and rare events.
<figure><img src="image_13.png" /><figcaption>Figure 12: Uncertainty evolution in time.</figcaption></figure>

Furthermore, as expected, higher uncertainty is observed in the transient part of the signal, compared to the limit cycle. This phenomenon is formally illustrated in Figure 12. To analyse temporal variations in uncertainty for a fixed parameter, we compute the spatially averaged uncertainty $\nu_t$ at each time step $t$ over $n_{xy}$ spatial locations,

$$ \nu_t = \frac{1}{n_{xy}} \sum_{d=1}^{n_{xy}} \nu_{d,t} \qquad (12) $$

where,

$$ \nu_{d,t} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (\phi_{i,d,t} - \bar{\phi}_{d,t})^2}. \qquad (13) $$

The width of the bell curve in Figure 12 aligns with the transient phase of the dynamic evolution visible on the kinetic energy signal. Starting from the initial condition, the uncertainty increases as the system transitions towards the limit cycle. Once this stable state is reached, the uncertainty decreases and stabilises at a plateau. This can be explained by two factors: (1) the regime transition induces variability, which like any other source of variability, is associated with lower confidence, and (2) the training distribution contains fewer transient states compared to the limit cycle. As a result, the model has been trained more extensively on the stable state than on the transient phases.

Thus, we highlight two key factors driving uncertainty: the system's inherent variability and the distance from the training distribution. While we have limited control over the first factor, we can influence the second, making it particularly valuable for UQ-driven adaptive training, as demonstrated in Section 5.

### 4.3. Uncertainty in the Parameter space

Uncertainty Quantification also provides a confidence measure across the parameter space, aggregating the spatial and temporal components,

$$ \nu_{\xi} = \frac{1}{n_{xy} n_t} \sum_{d=1}^{n_{xy}} \sum_{t=1}^{n_t} \nu_{\xi, d, t}, \qquad (14) $$

where,

$$ \nu_{\xi, d, t} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (\phi_{i, \xi, d, t} - \bar{\phi}_{\xi, d, t})^2}, \qquad (15) $$
with $n_t$ the time discretisation and $\bar{\phi}_{\epsilon,d,t}$ the mean of the ensemble of size $n$.

As shown in Figure 8, the model can learn and infer within the 2D parameter space. Specifically, it can interpolate within the training subspace and extrapolate with varying accuracy depending on the region of the parameter space. The further the inference distribution deviates from the training distribution, the greater the uncertainty of the model, although not uniformly. The direction in which the inference distribution moves away from the training distribution affects the performance of the model. Our uncertainty quantification method can measure this non-Euclidean distance. Figure 13 illustrates how our uncertainty quantification process allows us to identify these regions during inference, thus reducing the need for simulation data for evaluation. To improve the performance of the model, it could be retrained using samples from and only from these uncertain regions. An important observation is the strong agreement between the a-posteriori performance evaluation shown in Figure 8 and the a-posteriori inference of these performances. Both highlight the same at-risk dynamic regime. Therefore, our UQ framework is useful for fine-tuning, adaptive training, or error detection.

<figure><img src="image_14.png" /><figcaption>Figure 13: Scaled uncertainty in the $(R_e, \beta)$ parameter space.</figcaption></figure>

In particular, the model exhibits lower uncertainty near the training samples, further confirming the correlation between confidence and distance from the training distribution. Furthermore, the model shows particularly low confidence near the bifurcation boundary, demonstrating its strong dynamic awareness. Even without training in the transition regime, the model detects this transition and recognises it as a potential challenge.

We can use this parameter-dependent measure of uncertainty to guide the selection of parameter sets for training. This approach supports adaptive sampling, allowing optimised learning and generalisation for a range of dynamic behaviours, while minimising the amount of training data required. In Section 5, we demonstrate a clear correlation between regions of poor model performance and high uncertainty, and identify optimal areas for retraining.
### 4.4. Robustness across the parametric regime

The performance of the probabilistic model can be assessed using the Relative Mean Squared Error (MSE) calculated using test data, as

$$ \text{Relative MSE}_{\xi} = \frac{\sum_{d,t} (\hat{\Phi}_{\xi,t,d} - \Phi_{\xi,t,d})^2}{\sum_{d,t} \Phi_{\xi,d,t}^2} \times 100\% $$

where, $\hat{\Phi}_{\xi,d,t}$ is the predicted value for parameter $\xi$ at time $t$ and dimension $d$ and $\Phi_{\xi,d,t}$ is the ground truth value for parameter $\xi$ at time $t$ and spatial location $d$. Another measure is the Continuous Rank Probability Score (CRPS) computed for an ensemble of $n$ sampled trajectories given by:

$$ \text{CRPS}_{\xi} = \frac{1}{n} \sum_{i}^{n} (\hat{\Phi}_{i,\xi,d,t} - \Phi_{\xi,d,t})^2 - \frac{1}{2n^2} \sum_{i}^{n} \sum_{j \neq i}^{n} (\hat{\Phi}_{i,\xi,d,t} - \hat{\Phi}_{j,\xi,d,t})^2. $$

For a model parametrised by the Reynolds number ($R_e$) and trained for $R_e = 90$, $R_e = 120$, these two measures are reported in table 1.

<table><caption>Table 1: Performance metrics evaluated on test data across various Reynolds numbers.</caption><tr><td>R<sub>e</sub> :</td><td>50</td><td>60</td><td>70</td><td>80</td><td><b>90</b></td><td>100</td><td>110</td><td><b>120</b></td><td>130</td><td>140</td></tr><tr><td>Relative MSE (%)</td><td>0.24</td><td>0.20</td><td>1.56</td><td>1.61</td><td><b>0.28</b></td><td>0.71</td><td>2.49</td><td><b>0.42</b></td><td>7.57</td><td>9.09</td></tr><tr><td>CRPS (×10<sup>-5</sup>)</td><td>6.3</td><td>5.0</td><td>7.17</td><td>6.65</td><td><b>5.37</b></td><td>3.51</td><td>6.89</td><td><b>3.22</b></td><td>13.7</td><td>24.5</td></tr></table>

The table indicates that with just two training points, the relative error stays below 10%, and drops even further to below 2.5% when the Reynolds number $R_e$ is less than 130, which is within the range of the training distribution. This trend persists even when extrapolating outside the training distribution, particularly for $R_e < 90$. As a result, the model maintains satisfactory predictive performance, even in critical applications, highlighting its robustness. Additionally, these findings suggest that the impact of extrapolation on performance is not uniform across the entire parameter space.

The Continuous Ranked Probability Score (CRPS) is on the order of $10^{-4}$, indicating a very low value. Conceptually, CRPS measures the difference between the mean distance of the ensemble from the true target and the variance within the ensemble. A low CRPS indicates that each sampled trajectory is close to the true target trajectory and that the ensemble members themselves are close to each other. Since the model performs well near the training distribution, there is little variation in the ensemble predictions. However, as performance decreases (e.g., when $R_e > 130$), the CRPS increases, reflecting a greater variety within the ensemble. Such a low CRPS is synonymous with an accurate and reliable model.

the performance of a model that has been parametrised by the Reynolds Number ($R_e$) and trained for $R_e = 50$, $R_e = 60$, $R_e = 70$, $R_e = 80$, $R_e = 90$, $R_e = 120$, $R_e = 140$ is given in table 2, which illustrates the impact of retraining on the model's performance across the parameter space.

<table><caption>Table 2: Performance metrics evaluated on test data across various Reynolds numbers, with training data encompassing a broader range of Reynolds numbers.</caption><tr><td>R<sub>e</sub> :</td><td><b>50</b></td><td><b>60</b></td><td><b>70</b></td><td><b>80</b></td><td><b>90</b></td><td>100</td><td>110</td><td><b>120</b></td><td>130</td><td><b>140</b></td></tr><tr><td>Relative MSE (10<sup>-2</sup>%)</td><td><b>0.02</b></td><td><b>0.13</b></td><td><b>2.37</b></td><td><b>3.59</b></td><td><b>3.58</b></td><td>5.99</td><td>2.74</td><td><b>5.40</b></td><td>11.88</td><td><b>18.24</b></td></tr><tr><td>CRPS (×10<sup>-5</sup>)</td><td><b>0.32</b></td><td><b>1.09</b></td><td><b>2.61</b></td><td><b>5.59</b></td><td><b>3.32</b></td><td>3.04</td><td>4.87</td><td><b>5.98</b></td><td>11.01</td><td><b>13.23</b></td></tr></table>

The results show that expanding the training distribution reduces the relative MSE to below 0.2% on test data for the entire parameter range of interest, which, to our knowledge, is among the lowest relative MSE reported for the similar benchmarks of the “flow around a cylinder”. The CRPS remains low, but has not decreased as drastically as the relative MSE, indicating that while the model has become more accurate, the diversity of the ensemble has increased. As the dynamic range in the training distribution expands, the
model becomes aware of a wider range of possible dynamics, leading to a decrease in overall confidence. However, in Section 5 we show that the aggregated uncertainty is not as relevant as its distribution over the parameter space.

## 5. Adaptive Sampling for Optimised Generalisation

A parametrised model naturally prompts the question of how to select the appropriate parameters to ensure effective learning and inference across the entire parameter space. Since gathering training data for the full domain is often infeasible, we are typically limited to a small subset of the parameter space. The central issue, therefore, is identifying the minimal set of parameters that will offer the most comprehensive information, enabling the model to generalise effectively across the entire range of parameters. We aim to use Uncertainty Quantification (UQ) to efficiently explore the parameter space and construct an adaptive sampling strategy. To this end, the high uncertainties identify regions where model performance is suboptimal and indicate the region of the parameter space where data need to be collected. This approach leads to an adaptive sampling technique that adjusts data sampling based on model performance or uncertainty to improve learning efficiency and accuracy. Our approach is implemented during inference of the UQ framework and reduces the reliance on prior knowledge of the parameter space and on the need for validation data throughout the adaptive procedure, which would be necessary with performance/error driven adaptive sampling.

This approach can only be effective if the uncertainty can serve as a reliable measure of model performance or predictability. To validate this, we first examine whether there is a strong correlation between the *a posteriori* error and the *a priori* uncertainty. To do this, we choose a parametrised model with a variable Reynolds number $\xi = Re$ and a fixed aspect ratio ($\beta = 1.35$). This configuration leads to a bifurcation at Reynolds number $\approx 60$. The dynamical regime before this critical Reynolds number ($Re_c$) shows a stationary behaviour, while the regime after the bifurcation results in a time-varying limit cycle, as illustrated in Figure 4 and in Figure 14 by the different colours of the shaded regions. The experiment is performed over a parameter range from $Re = 50$ in the pre-bifurcation regime to $Re = 140$ in the post-bifurcation regime. The model infers over this entire range on 10 equally spaced Reynolds number values, for which validation data are available. We obtain 10 measurements of $\nu_\xi$ that we compare with the model's performance over the same range. Performance is evaluated using the Mean Squared Error (MSE), scaled by the order of magnitude of the signal over the range. Note that in real-world scenarios, validation data over the full range would not be necessary. The scaled MSE is computed here *a posteriori* solely to assess and validate our uncertainty quantification.

At initialisation, the model is trained with only two samples in the post-bifurcation regime: $Re = 90$ and $Re = 120$, exposing the model to minimal parameter variety. Figure 14-a shows the initial uncertainty distribution over the parameter range $\nu_\xi$, while Figure 14-b shows the scaled error. There is a strong correlation between the two values as formally shown in Figure 14-c, with the linear fit between $\nu_\xi$ and the error, along with the corresponding Pearson correlation index. In this first step, we observe a Pearson index of $r = 0.87 > 0.5$, indicating a strong, positive linear correlation. Not surprisingly, the model is most confident for the in-distribution parameters and specifically for the two training samples. It is least confident in the pre-bifurcation regime. Since the model has never encountered this regime, it exhibits both greater uncertainty and poorer performance. This analysis confirms that the UQ measure serves as a reliable indicator of model predictability or error, and provides the basis for an adaptive sampling strategy that iteratively guides the retraining of the model until uncertainty is minimised over the parameter range of interest.

The uncertainty distribution given by $\nu_\xi$ in Figure 14-a shows the highest uncertainty at $Re = 50$. Therefore, we feed the model with this data and perform the first retraining. This retraining is significantly shorter than the initial training, with fewer initial samples included to avoid forgetting. After retraining, both uncertainty and MSE are recomputed at inference and *a posteriori*, respectively, and exhibit a complete redistribution across the parameter range, as shown in Figure 14-d,e.

After the first retraining, the confidence in the pre-bifurcation segment improves significantly. The stationary-like dynamics make inference easier for the model as the output does not change with time. At
this stage, error and uncertainty remain well correlated, as shown in Figure 14-f. The highest uncertainty now
shifts to the other side of the bifurcation, at R<sub>e</sub> = 80, suggesting the need to refine the model's understanding
of the segment near the critical Reynolds number. The data associated with this parameter value becomes
the next candidate point for retraining. Interestingly, at this stage, extrapolation beyond the training
distribution (R<sub>e</sub> > 120) is less challenging in terms of both uncertainty and error compared to the lower range
(R<sub>e</sub> < 90). The retraining process can be extended by iteratively adding samples with the highest uncertainty
until the uncertainty converges below a desired value. The complete iterative retraining process, together
with the uncertainty distribution ν<sub>ξ</sub>, the scaled MSE and their correlations is shown in Figure 14. The results
throughout this process indicate that points with the highest Mean Squared Error (MSE) consistently align
with those exhibiting the highest uncertainty. A clear and consistent linear relationship between error
and uncertainty is evident. Following each retraining iteration, the Pearson correlation coefficient between
the UQ and MSE exceeds 0.7, underscoring a strong positive linear correlation. This demonstrates the
effectiveness of our uncertainty quantification method in anticipating poor performance across the parameter
space and for adaptive sampling.
<figure><img src="image_15.png" /><figcaption>Iteration 5: Retraining at Reynolds number 60</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14: Uncertainty, scaled MSE and their correlation throughout a retraining process</figcaption></figure>
After completing five retraining iterations, the UP-dROM method demonstrates consistent and uniform learning performance across the entire range of model parameters, effectively capturing the underlying bifurcation dynamics. Notably, the model’s uncertainty appears to have converged to a stable final distribution, denoted as $v_ξ$. The model progressively identifies parameter points that minimise the error at each training step. As retraining continues at these optimal points, we observe a convergence of both the scaled mean squared error (MSE) and the uncertainty, which both tend towards an optimal distribution.

# 6. Conclusions and perspectives

This paper presents a robust and parametrised reduced-order model (ROM) for time-dependent partial differential equations applied to the Navier-Stokes equations in the laminar regime. The model uses memory embedding and variational inference, first, a nonlinear variational autoencoder is used to reduce the high-dimensional dynamics to a low-dimensional latent space and second, an attention-based transformer captures time dependencies in the reduced manifold. In addition, cross-attention incorporates dependencies on external parameters. This allows UP-dROM to handle different dynamical regimes, making it adaptable to external excitation variables and enhancing generalisation.

Our model demonstrates high accuracy, with a relative mean square error consistently below 0.5% on the test set for seen parameters. Interpolation between seen parameters remains below 2.5%, even with limited parameter diversity. Extrapolation is effective in some regions, but less reliable in others. However, the model’s Uncertainty Quantification (UQ) framework allows these uncertain regions to be identified a priori. A key contribution of this work is the integration of a computationally efficient UQ framework that provides instantaneous confidence estimates without adding to the computational cost during inference. The approach links uncertainty to both system variability and distance from the training distribution, providing an interpretable measure of confidence across time, spatial domain, and parameter space.

We show that UQ is not only valuable in itself, but also a useful tool for efficiently sampling the parameter space. Since collecting training data across the entire domain is often impractical, training is restricted to a small subset of parameters. Furthermore, we show a strong positive linear correlation between the error and the uncertainty in the parameter space. Therefore, quantifying the uncertainty of the model can guide adaptive sampling, focusing retraining efforts on regions of high uncertainty. This reduces the need for validation data across the entire domain, optimising learning and enhancing generalisation across different dynamical regimes. In addition, the framework effectively captures bifurcation dynamics and complex parametric dependencies, making it well suited for modelling nonlinear fluid systems.

Future work will explore data assimilation strategies to refine the model in uncertain regions. The model can identify regions of high uncertainty in physical space where sparse sensor measurements can enhance the model without requiring additional full-state training data. In addition, using the model’s parametrisation as a control variable within a control framework could open new possibilities for flow optimisation.
References

[1] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational physics 378 (2019) 686-707.

[2] S. Cai, Z. Mao, Z. Wang, M. Yin, G. E. Karniadakis, Physics-informed neural networks (pinns) for fluid mechanics: A review, Acta Mechanica Sinica 37 (12) (2021) 1727-1738.

[3] M. Cranmer, S. Greydanus, S. Hoyer, P. Battaglia, D. Spergel, S. Ho, Lagrangian neural networks (03 2020). doi:10.48550/arXiv.2003.04630.

[4] M. Mattheakis, D. Sondak, A. S. Dogra, P. Protopapas, Hamiltonian neural networks for solving equations of motion, Phys. Rev. E 105 (2022) 065305. doi:10.1103/PhysRevE.105.065305.
URL https://link.aps.org/doi/10.1103/PhysRevE.105.065305

[5] R. T. Chen, Y. Rubanova, J. Bettencourt, D. K. Duvenaud, Neural ordinary differential equations, Advances in neural information processing systems 31 (2018).

[6] K. He, X. Zhang, S. Ren, J. Sun, Deep residual learning for image recognition, 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (2015) 770–778.
URL https://api.semanticscholar.org/CorpusID:206594692

[7] N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, A. Anandkumar, Neural operator: Learning maps between function spaces with applications to pdes, Journal of Machine Learning Research 24 (89) (2023) 1–97.

[8] K. Azizzadenesheli, N. Kovachki, Z. Li, M. Liu-Schiaffini, J. Kossaifi, A. Anandkumar, Neural operators for accelerating scientific simulations and design, Nature Reviews Physics 6 (5) (2024) 320–328.

[9] B. Bahmani, S. Goswami, I. G. Kevrekidis, M. D. Shields, A resolution independent neural operator, arXiv preprint arXiv:2407.13010 (2024).

[10] V. Shinde, D. Gaitonde, Galerkin-pod reduced-order modeling for perturbation analysis and sparse state estimation of compressible flows, 2021. doi:10.2514/6.2021-2603.

[11] R. Pinnau, Model Reduction via Proper Orthogonal Decomposition, Springer Berlin Heidelberg, Berlin, Heidelberg, 2008, pp. 95–109. doi:10.1007/978-3-540-78841-6_5.
URL https://doi.org/10.1007/978-3-540-78841-6_5

[12] P. Schmid, J. Sesterhenn, Dynamic mode decomposition of numerical and experimental data, Journal of Fluid Mechanics 656 (11 2008). doi:10.1017/S0022112010001217.

[13] P. J. SCHMID, Dynamic mode decomposition of numerical and experimental data, Journal of Fluid Mechanics 656 (2010) 5–28. doi:10.1017/S0022112010001217.

[14] P. J. Baddoo, B. Herrmann, B. J. McKeon, J. Nathan Kutz, S. L. Brunton, Physics-informed dynamic mode decomposition, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 479 (2277) (2023) 20220576. doi:10.1098/rspa.2022.0576.
URL http://doi.org/10.1098/rspa.2022.0576

[15] G. Haller, S. Ponsioen, Nonlinear normal modes and spectral submanifolds: existence, uniqueness and use in model reduction, Nonlinear Dynamics 86 (3) (2016) 1493–1534. doi:10.1007/s11071-016-2974-z.
URL https://doi.org/10.1007/s11071-016-2974-z
[16] G. Buza, S. Jain, G. Haller, Using spectral submanifolds for optimal mode selection in nonlinear model reduction, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 477 (2021) 20200725. doi:10.1098/rspa.2020.0725.
URL http://doi.org/10.1098/rspa.2020.0725

[17] K. Lee, K. T. Carlberg, Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders, J. Comput. Phys. 404 (C) (Mar. 2020). doi:10.1016/j.jcp.2019.108973.
URL https://doi.org/10.1016/j.jcp.2019.108973

[18] Özalp Elise, M. Luca, Stability analysis of chaotic systems in latent spaces, Nonlinear Dynamics (2025). doi:10.1007/s11071-024-10712-w.
URL https://doi.org/10.1007/s11071-024-10712-w

[19] E. Menier, M. A. Bucci, M. Yagoubi, L. Mathelin, M. Schoenauer, Cd-rom: Complemented deep - reduced order model, Computer Methods in Applied Mechanics and Engineering 410 (2023) 115985. doi:https://doi.org/10.1016/j.cma.2023.115985.
URL https://www.sciencedirect.com/science/article/pii/S0045782523001081

[20] S. L. Brunton, M. Budišić, E. Kaiser, J. N. Kutz, Modern koopman theory for dynamical systems, SIAM Review 64 (2) (2022) 229-340. arXiv:https://doi.org/10.1137/21M1401243, doi:10.1137/21M1401243.
URL https://doi.org/10.1137/21M1401243

[21] E. Menier, S. Kaltenbach, M. Yagoubi, M. Schoenauer, P. Koumoutsakos, Interpretable learning of effective dynamics for multiscale systems (Submitted).
URL https://arxiv.org/abs/2309.05812

[22] P. Gupta, P. J. Schmid, D. Sipp, T. Sayadi, G. Rigas, Mori-zwanzig latent space koopman closure for nonlinear autoencoder, ArXiv abs/2310.10745 (2023).
URL https://api.semanticscholar.org/CorpusID:264172182

[23] A. Kristiadi, M. Hein, P. Hennig, Learnable uncertainty under laplace approximations, in: C. de Campos, M. H. Maathuis (Eds.), Proceedings of the Thirty-Seventh Conference on Uncertainty in Artificial Intelligence, Vol. 161 of Proceedings of Machine Learning Research, PMLR, 2021, pp. 344-353.
URL https://proceedings.mlr.press/v161/kristiadi21a.html

[24] T. Simpson, K. Vlachas, A. Garland, N. Dervilis, E. Chatzi, Vprom: a novel variational autoencoder-boosted reduced order model for the treatment of parametric dependencies in nonlinear systems, Scientific Reports 14 (03 2024). doi:10.1038/s41598-024-56118-x.

[25] A. Solera-Rico, C. Sanmiguel Vila, M. Gómez-López, Y. Wang, A. Almashjary, S. T. M. Dawson, R. Vinuesa, β-variational autoencoders and transformers for reduced-order modelling of fluid flows, Nature Communications 15 (1) (2024) 1361. doi:10.1038/s41467-024-45578-4.
URL https://doi.org/10.1038/s41467-024-45578-4

[26] A. Damianou, Deep gaussian processes and variational propagation of uncertainty, Ph.D. thesis, University of Sheffield (2015).

[27] T. G. Dietterich, Ensemble methods in machine learning, in: Multiple Classifier Systems, Springer Berlin Heidelberg, Berlin, Heidelberg, 2000, pp. 1-15.

[28] Y. Gal, Z. Ghahramani, Dropout as a bayesian approximation: Representing model uncertainty in deep learning, in: M. F. Balcan, K. Q. Weinberger (Eds.), Proceedings of The 33rd International Conference on Machine Learning, Vol. 48 of Proceedings of Machine Learning Research, PMLR, New York, New York, USA, 2016, pp. 1050-1059.
URL https://proceedings.mlr.press/v48/gal16.html
[29] I. Shokar, R. Kerswell, P. Haynes, Stochastic latent transformer: Efficient modeling of stochastically forced zonal jets, Journal of Advances in Modeling Earth Systems 16 (06 2024). doi: 10.1029/2023MS004177.

[30] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. u. Kaiser, I. Polosukhin, Attention is all you need 30 (2017). URL https://proceedings.neurips.cc/paper_files/paper/2017/file/3f5ee243547dee91fbd053c1c4a845aa-Paper.pdf

[31] L. Pinheiro Cinelli, M. Araújo Marins, E. A. Barros da Silva, S. Lima Netto, Variational autoencoder, in: Variational methods for machine learning with applications to deep networks, Springer, 2021, pp. 111–149.

[32] S. Kneer, T. Sayadi, D. Sipp, P. Schmid, G. Rigas, Symmetry-aware autoencoders: s-pca and s-nlPCA, arXiv preprint arXiv:2111.02893 (2021).

[33] C. Schering, G. Rigas, D. Sipp, P. Schmid, T. Sayadi, An adaptive learning strategy for surrogate modeling of high-dimensional functions - application to unsteady hypersonic flows in chemical nonequilibrium, Computer Physics Communications 307, publisher Copyright: © 2024 Elsevier B.V. (Feb. 2025). doi:10.1016/j.cpc.2024.109404.

[34] B. Moore, Principal component analysis in linear systems: Controllability, observability, and model reduction, IEEE Transactions on Automatic Control 26 (1) (1981) 17–32. doi:10.1109/TAC.1981.1102568.

[35] A. Mavi, A. C. Bekar, E. Haghighat, E. Madenci, An unsupervised latent/output physics-informed convolutional-lstm network for solving partial differential equations using peridynamic differential operator, Computer Methods in Applied Mechanics and Engineering 407 (2023) 115944. doi:https://doi.org/10.1016/j.cma.2023.115944.
URL https://www.sciencedirect.com/science/article/pii/S0045782523000671

[36] Y. T. Lin, Y. Tian, D. Livescu, M. Anghel, Data-driven learning for the mori-zwanzig formalism: A generalization of the koopman learning framework, SIAM Journal on Applied Dynamical Systems 20 (4) (2021) 2558–2601. arXiv:https://doi.org/10.1137/21M1401759, doi:10.1137/21M1401759.
URL https://doi.org/10.1137/21M1401759

[37] F. Takens, Detecting strange attractors in turbulence, 1981.
URL https://api.semanticscholar.org/CorpusID:10823331

[38] I. Zighed, p-dynsys_encodertransformerdecoder, https://github.com/IsmaelZig-SU/p-DynSysENCoderTransformerDecoder (2025).

[39] M. F. de Pando, miguelfp/ibmos: Initial release (Apr. 2020). doi:10.5281/zenodo.3757783.
URL https://doi.org/10.5281/zenodo.3757783

[40] I. Zighed, Training and validation data of updrom - flow passed a 2d blunt obstacle, https://drive.google.com/file/d/1iZpAsPFqziRx3hTSnsfMqt9czlftcq3i/view?usp=sharing (2025).

[41] B. Jia, H. Wu, K. Guo, Chaos theory meets deep learning: A new approach to time series forecasting, Expert Systems with Applications 255 (2024) 124533. doi:https://doi.org/10.1016/j.eswa.2024.124533.
URL https://www.sciencedirect.com/science/article/pii/S0957417424014003

[42] N. Geneva, N. Zabaras, Transformers for modeling physical systems, Neural Networks 146 (2022) 272–289. doi:https://doi.org/10.1016/j.neunet.2021.11.022.
URL https://www.sciencedirect.com/science/article/pii/S0893608021004500
[43] P. J. Strykowski, K. R. Sreenivasan, On the formation and suppression of vortex ‘shedding’ at low reynolds numbers, Journal of Fluid Mechanics 218 (1990) 71–107. doi:10.1017/S0022112090000933.

[44] D. Canuto, K. Taira, Two-dimensional compressible viscous flow around a circular cylinder, Journal of Fluid Mechanics 785 (2015) 349–371. doi:10.1017/jfm.2015.635.
# 7. Appendix

## 7.1. Generalisation to Different Physics

We believe that UP-dROM is physics agnostic and could therefore apply to other family of equations, especially the The Kuramoto-Sivashinsky (KS) equation. The KS equation is a nonlinear partial differential equation used as a benchmark in turbulence modelling and the study of chaotic systems. It is given by:

$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \frac{\partial^2 u}{\partial x^2} + \nu \frac{\partial^4 u}{\partial x^4} = 0 $$

where $u(x, t)$ represents the field of interest, $x$ is the spatial coordinate, $t$ is time, and $\nu$ is a positive parameter. The KS equation is valued for its ability to capture the essential features of turbulence, such as chaotic behaviour and complex spatiotemporal patterns, while remaining relatively simple and analytically tractable, making it an ideal testbed for developing and validating UP-dROM through turbulent or near turbulent regimes. We compare the inferred and ground truth domains for different values of parameter $\nu$ in Figure 15
<figure><img src="image_17.png" /><figcaption>Figure 15: Kuramoto–Sivashinsky model performance comparing the ground truth $u$ and the prediction $\hat{u}$ across a range of $v$ values</figcaption></figure>

## 7.2. Model training and architecture

Transformer and autoencoder are trained together to promote a symbiotic relationship between the two elements, thus fostering an organic modal representation. The total loss is calculated considering a prediction horizon $h$ and a lookback window $q$. Tables 4 and 5 provide further details on the selection of hyperparameters across different test cases, which were chosen through grid search to balance computational efficiency and model accuracy. The loss function is defined as:

$$ \mathcal{L} = \frac{\lambda}{m} \left( \sum_{k=0}^{m} \|\phi_{n-q:n} - \mathcal{D}\mathcal{E}(\phi_{n-q:n})\| + \beta \cdot \text{KLD} \right) + \frac{1}{m} \sum_{k=0}^{m} \|z_{n+1:n+h} - \hat{z}_{n+1:n+h}\| + \frac{1}{m} \sum_{k=0}^{m} \|\phi_{n+1:n+h} - \mathcal{D}(\hat{z}_{n+1:n+h})\|. \quad (16) $$
where $\mathcal{E}$ is the encoding process, $\mathcal{D}$ the decoding process, KLD the Kullback-Leibler divergence and $\hat{z}$ the latent variables. In addition, $\lambda$ and $\beta$ are regularisation coefficients. The first term of this loss promotes a mirroring process between encoder and decoder, similar to a classical projection. The Kullback-Leibler divergence term is related to the use of the VAE and is explained in Section 4.

<table><caption>Table 3: Loss function parameters</caption><thead><tr><th></th><th>λ</th><th>β</th></tr></thead><tbody><tr><td>Navier Stokes</td><td>100</td><td>1e-4</td></tr><tr><td>KS</td><td>100</td><td>1e-4</td></tr></tbody></table>

$$ \mathrm{KLD} = -\frac{1}{2} \sum_{i=1}^{Z} (\sigma_i^2 + \mu_i^2 - 1 - \log(\sigma_i^2)) \qquad (17) $$

<table><caption>Table 4: Model Hyperparameters for Navier Stokes and KS</caption><thead><tr><th></th><th>Prediction Horizon</th><th>Lookback Window</th><th>Space (dim)</th><th>Time (dim)</th><th>Latent Dimensions</th></tr></thead><tbody><tr><td>Navier Stokes</td><td>10</td><td>10</td><td>[131 × 100]</td><td>3033</td><td>4</td></tr><tr><td>KS</td><td>30</td><td>30</td><td>255</td><td>3033</td><td>32</td></tr></tbody></table>

<table><caption>Table 5: Transformer Hyperparameters for Navier Stokes and KS</caption><thead><tr><th></th><th>Hidden Dimensions</th><th>Attention Heads</th><th>Attention Blocks</th></tr></thead><tbody><tr><td>Navier Stokes</td><td>64</td><td>8</td><td>1</td></tr><tr><td>KS</td><td>128</td><td>16</td><td>2</td></tr></tbody></table>