INTRODUCTION

Turbulence remains a central challenge in classical physics [1]. Its presence is ubiquitous in natural and engineered flows from sub-meter to planetary scales. Accurate prediction of turbulent flows has far-reaching impacts on modern society, ranging from reducing uncertainties in climate projections and improving fidelity in extreme weather events [2] to designing safer and more efficient mission-critical systems such as aircraft engines and power plants [3]. The chaotic and multi-scale nature of turbulence [4] is responsible for critical physics in these systems, such as momentum and heat transport in planetary boundaries and flow through engine turbines and pollutant transport in rivers and oceans. Despite the importance of turbulent transport in such processes, our current modeling of such effects in state-of-the-art simulations is still crude. For example, in weather forecast models, the vertical turbulent mixing is still represented as one-dimensional parameterization schemes known as the planetary boundary layer scheme [5]. Similarly, in industrial computational fluid dynamics, Reynolds-averaged Navier–Stokes (RANS) solvers are still the workhorse tool for routine simulations [6], where the RANS equations govern the mean flow field and the “averaged effects” of the turbulent fluctuations, introduced through averaging as the Reynolds stress, are represented with turbulence models [7,8]. This is because such applications require either long-time integration or many simulations for design optimization or uncertainty quantification. It becomes prohibitively expensive to fully resolve all the turbulence scales (as in direct numerical simulations) or even to partially resolve them (as in large eddy simulations). Despite the growth of computational resources, flow solvers based on pure modeling or parameterization of the turbulence (referred to as turbulence modeling collectively hereafter) will remain as the backbone in many fields in the years to come. Therefore, the theoretical foundation and practical development of turbulence will continue to have societal and scientific importance.

Pursuit of a universal constitutive relation

More than half a century ago, Lumley [9] postulated the existence of a turbulent constitutive relation that maps the mean velocities to the Reynolds stresses. From a phenomenological perspective, such a relation shall be inferred from data and constrained by mathematical invariance and symmetry, which in principle enables transferability across similar flows. Mathematically, it defines the Reynolds stress at a given location as a functional of the mean velocity field evaluated over a surrounding neighborhood, which may in the limiting case reduce to a purely local dependence. Using experimental data, Lumley established that such a local constitutive relation exists for homogeneous shear flows and homogeneous strain, both of which behave analogously to viscoelastic fluids with an increasing relaxation timescale [9]. Despite its conceptual depth, this hypothesis has had limited impact on the subsequent development of turbulence modeling. Although numerous turbulence models have since been proposed, none has achieved the universality across a broad range of flows envisioned by Lumley. Recent analyses have further emphasized fundamental limitations of single-point, local closures in representing non-local transport and non-equilibrium effects in turbulence [10], as purely local constitutive relations cannot capture the history dependence and spatial interactions inherent in turbulent flows. Consequently, the existence of a universal, local constitutive relation remains uncertain under practical modeling constraints.

From universality to unification

Rather than pursuing a universal constitutive relation of uncertain existence, this work focuses on developing a unified turbulence model that can represent multiple flow mechanisms within a single framework. Here, a unified model denotes a closure capable of handling multiple flow mechanisms, such as attached boundary layers, separated flows, and secondary flows, without manual zoning, switching, or blending [11,12]. Such unification is particularly attractive for industrial applications such as aircraft aerodynamics and turbine flows, where multiple flow mechanisms coexist, interact, and cannot be clearly separated spatially. Developing a unified foundation model of this type is the central objective of this work. Nevertheless, even a more limited degree of unification, referred to here as a specialist model, would already offer substantial practical value [13,14].

Efforts toward a unified turbulence model

Achieving a unified turbulence model relies on training with diverse flows that capture multiple flow mechanisms. A number of strategies have been pursued with varying degrees of success. For example, multi-case training using symbolic model representations has produced models that perform well across several wall-bounded flows, representing some of the earliest steps toward unified turbulence modeling [15,16]. Nevertheless, the restrictive functional forms and reliance on derivative-free optimization can limit scalability, particularly when many competing objectives arise. An alternative strategy, which we refer to as “unification through aggregation”, relies on expert models tailored to individual flow mechanisms, with a classifier assigning one or more models to each spatial location during prediction [17,18]. Such approaches, however, require multiple model evaluations and can struggle when different flow mechanisms strongly interact within the same region. Yet, another strategy toward a unified model is to provide users with a generalized model that can be tuned for individual flows, exemplified by the generalized k–ω model (GEKO) [19]. This model relies on manual tuning for different flows, which limits its ability to achieve true unification. Similar efforts toward a generalizable closure have also emerged in wall-modeled large eddy simulations, such as the knowledge-integrated additive learning [20] and building block approaches [21]. However, to this day, a unified RANS turbulence model that can seamlessly handle coexisting, vastly different flow mechanisms in a single flow (for example, flows with massive separation and swirling simultaneously) does not yet exist. The difficulty lies in the fact that different flow regimes impose conflicting requirements on the model: adjustments that improve one regime can degrade another. This conflict fundamentally limits the transferability of existing models across diverse flow regimes. As noted in the NASA Langley Turbulence Modeling Symposium [11], after a decade of intense research, the community has yet to produce a data-driven turbulence model that exceeds current models in terms of predictability, generality, and robustness.

Unification via multi-objective learning

In this work, we develop a unified data-driven turbulence modeling framework that learns a single closure capable of representing multiple flow mechanisms. The unification is achieved through a strategy that balances conflicting objectives and branches across non-conflicting ones. As illustrated in Fig. 1, the proposed machine-learning framework proceeds through three steps that together yield a unified turbulence model. First, we employ a physically consistent, frame-invariant model representation that enables a single turbulence model to adapt to different flow mechanisms without manual switching (Fig. 1a) [22,23]. Second, we use a distance-based selection strategy to automatically identify a compact and representative set of training flows that spans the relevant flow physics while avoiding redundancy (Fig. 1b). Finally, we devise an ensemble-based, multi-objective learning strategy to balance competing objectives across flow regimes and, through this process, learn a unified turbulence model from heterogeneous flows and sparse, indirect observations [24] (Fig. 1c). The resulting foundation model captures multiple flow mechanisms within a single set of network weights and is seamlessly integrated into a Reynolds-averaged Navier–Stokes solver in OpenFOAM.

Figure 1.: The proposed framework proceeds through three steps to construct a unified turbulence model. The learning is formulated as a multi-objective optimization problem, yielding a single neural-network-based model that reconciles competing objectives across flows and quantities of interest. (a) First, a physically consistent and frame-invariant model representation is employed, in which the turbulent constitutive relation and transport equations are learned in a coupled and internally consistent manner under physical constraints, enabling a single model to adapt to different flow mechanisms without manual switching. (b) Second, a comprehensive dataset of flows is compiled, and a distance-based training-set selection strategy is used to automatically identify a compact and representative set of training cases by comparing probability distributions of local, frame-invariant flow features, thereby spanning relevant flow physics while avoiding redundancy. (c) Finally, an ensemble-based, multi-objective learning framework is applied to learn a unified model from diverse flows and sparse, indirect observations, balancing competing objectives across flows and quantities of interest. Taken together, these three components yield a unified foundation turbulence model that captures multiple flow regimes within a single set of learned network weights. For application-specific accuracy, the unified foundation model can be adapted into a specialist model through additive fine-tuning (see Fig. S1 in Supplementary Material).

Building on prior work that enables learning a turbulence model from indirect data [24], this study introduces several advances that together enable robust unification across flow regimes. First, we develop a flexible and physically consistent turbulence model representation in which the turbulent constitutive relation and transport equations are learned in a coupled manner, ensuring coherent closure behavior and turbulence scales [23]. Second, we introduce an automatic, distribution-based strategy for selecting representative training flows, coupled with the compilation of a comprehensive library of 36 canonical-to-complex flows—the most extensive dataset used to train a single turbulence model to date. Finally, and most importantly, we formulate turbulence model learning as a multi-objective optimization problem, allowing competing objectives arising from different flows and quantities of interest to be reconciled within a single unified model.

In particular, this formulation naturally accommodates multiple quantities of interest from the same flow (for example, drag and lift for an airfoil), which can introduce competing demands on the turbulence model (see Section S3.2.2).

Flexible and physically consistent model representation

We construct a turbulence model that combines high expressive power with physical consistency and numerical robustness by embedding invariance, structural hierarchy, and physics-based constraints directly into the model representation (see Fig. 1a). The tensor basis neural network (TBNN) framework [22] provides a flexible and frame-invariant realization of the general eddy-viscosity model, making it well suited for unified turbulence modeling across diverse flow regimes. In its standard form, the TBNN predicts coefficients for tensor bases of all polynomial orders through a single, shared network, implicitly weighting linear and higher-order nonlinear terms equally. This often introduces unnecessary nonlinearity that degrades generalization and numerical stability. To impose a clear physical structure, we introduce a parallel TBNN architecture (see Fig. S2b in Supplementary Material), in which shared invariant features are mapped to separate, order-specific network branches that predict coefficients for low-order and higher-order tensor bases independently, enabling physics-informed regularization of nonlinear contributions. Beyond the constitutive relation alone, turbulence modeling requires internal consistency with the turbulence transport equations that determine turbulent time and length scales. Rather than modifying either component in isolation [22], we treat coefficients in both the constitutive relation and the transport equations as learnable fields within a unified network architecture and optimize them in a coupled manner [23]. Finally, we impose physics-based constraints to ensure correct behavior in canonical limits, including decaying homogeneous isotropic turbulence, equilibrium homogeneous shear, compatibility with the logarithmic law of the wall, and edge of turbulence region [12,25]. Together, these design choices establish a turbulence model representation that is flexible yet structured, enforces internal consistency between constitutive and transport closures, and embeds physical constraints to ensure correct behavior in canonical flows. Despite these improvements, the model remains a single-point, local closure and thus inherits fundamental limitations. Such closures approximate inherently multi-point and non-local turbulence dynamics using local flow quantities and cannot distinguish between different flow structures that produce similar Reynolds stresses [26]. The present work aims to improve predictive accuracy within this framework rather than to overcome these intrinsic limitations.

Distribution-based training set selection

Developing a unified turbulence model requires training data that span diverse flow mechanisms while keeping the training cost manageable. To this end, we compiled a comprehensive library of benchmark flows from the literature, including both classical datasets [27] and recent datasets designed for data-driven turbulence modeling [28]. The resulting library comprises 36 flow cases spanning canonical to complex configurations; nine representative cases are selected for training using a hierarchical clustering strategy based on probability-distribution distances (described below), while the remaining cases are reserved for testing (see Fig. 2a). As extrapolation beyond the training feature space remains challenging for data-driven turbulence models [29,30], the evaluation focuses on generalization within the feature space spanned by the training flows. For each flow, appropriate quantities of interest (for example, velocity profiles, drag, and lift; see Table 1) are identified based on the underlying physics and modeling challenges, providing a rigorous basis for model evaluation. Unlike data-rich fields such as weather forecasting or genomics, turbulence modeling encompasses many distinct flow mechanisms but offers limited accessible data, making principled case selection essential and requiring an objective notion of similarity between flows. To define such distances systematically, we refrain from heuristic judgment, which lacks objectivity and reproducibility, and from pointwise flow comparisons, which are sensitive to mesh resolution. Instead, we assess flow similarity by comparing probability distributions of local, frame-invariant flow features (see Fig. 1b). Together, this dataset compilation and distribution-based selection strategy yields a training set that spans dominant flow mechanisms while preserving physical representativeness.

Figure 2.: Overview of the training and evaluation cases and performance of the unified foundation turbulence model. (a) Library of canonical and complex flows used for training and evaluation, comprising 36 cases in total. Canonical flows are grouped into four categories—attached boundary layers, free-shear flows, secondary flows, and separated flows—while three complex three-dimensional flows involving multiple mechanisms are shown in the inner circle. Nine representative cases (highlighted with filled circles) form the training set, and the remaining 27 cases are used for testing. (b) Radar chart comparing normalized misfits of the baseline and unified foundation models across flow categories. Smaller radial distances indicate improved performance. The unified foundation model demonstrates robust performance across all categories while maintaining accuracy comparable to the baseline for attached boundary layers.

Categories | Flows | Control parameter | Observations | No. Train | No. Test
Attached boundary layers | Flat plate |  | Velocities, friction |  | 1
 | Plane channel | Reynolds number | Velocities |  | 3
 | Airfoil (attached) | Shape, angle of attack | Lift, drag |  | 4
Separated flows | Curved step |  | Velocities | 1 | 
 | Hump |  | Velocities | 1 | 
 | Bump | Height | Friction | 1 | 4
 | Periodic hill | Slope steepness | Velocities | 1 | 3
 | Airfoil (separated) | Shape, angle of attack | Lift, drag | 1 | 3
Secondary flows | Square duct | Reynolds number | Velocities | 2 | 2
 | Rectangular duct | Aspect ratio | Velocities | 1 | 3
Free-shear flows | Round jet |  | Velocities | 1 | 
Complex 3D flows | Generic car |  | Drag |  | 1
 | 3D diffuser |  | Velocities, friction |  | 1
 | Generic aircraft | Angle of attack | Lift, drag, friction, pressure |  | 2
Totals | 9 | 27

Multi-objective ensemble learning from sparse, indirect observations

To learn a unified turbulence model from sparse, indirect observations across multiple flows, we develop a learning framework that combines ensemble-based inference with multi-objective optimization (see Fig. 1c). In practice, most experimental and engineering datasets provide only indirect measurements, such as velocity profiles or integral forces, rather than full-field turbulence quantities [24,31]. Learning from such data using gradient-based methods typically requires adjoint-enabled or fully differentiable solvers [32–34], which are often unavailable in complex engineering applications. We therefore employ a regularized ensemble Kalman learning framework [35], which enables non-intrusive parameter updates through ensemble-based covariances between model predictions and observations while constraining deviations from a physically meaningful baseline model [36]. Beyond learning from indirect data, unified turbulence modeling requires balancing performance across multiple flows and quantities of interest. Rather than optimizing a single aggregated loss, we formulate training as a multi-objective optimization problem and seek Pareto-optimal solutions that reconcile competing objectives across flow regimes [37,38]. Although originally developed for gradient-based learning, this framework naturally extends to ensemble-based updates by interpreting objective-specific ensemble corrections as generalized descent directions. Together, this ensemble-based, multi-objective formulation enables robust learning of a single turbulence model across diverse flows and observational constraints. Sensitivity analyses of the training set selection and model size further confirm this robustness (Section S3).

When application-specific accuracy is prioritized, the unified foundation model can be further adapted into a specialist model through additive fine-tuning, in which a compact correction module is trained while all foundation parameters remain fixed. This minimal-modification strategy reallocates model capacity toward targeted flow mechanisms while preserving the generalization and robustness inherited from the unified foundation model.

In summary, the proposed framework integrates a flexible and physically consistent model representation with representative training-case selection and multi-objective ensemble learning. The framework targets generalization within the feature space spanned by the training flows, while regularization toward the baseline model ensures stable behavior outside this regime. As demonstrated in the Results section, this integrated design enables robust and generalizable predictions across diverse flow phenomena without the need for manual zoning or handcrafted blending. Together, these elements advance turbulence modeling toward solutions that are both predictive and practical for large-scale engineering applications. More importantly, the methodology itself is scalable to a large number of flows and target quantities of interest—up to 40 objectives [37]—offering a viable path toward unifying diverse benchmark flow mechanisms and enabling application to truly complex industrial configurations.