# Data-driven discovery of dimensionless numbers and governing laws: A scaling technique for nuclear test facility design

Camille Razaire, Amandine Marrel, Bertrand Iooss, Sébastien Renaudière de Vaux, Jérome Cardolaccia

► To cite this version:

Camille Razaire, Amandine Marrel, Bertrand Iooss, Sébastien Renaudière de Vaux, Jérome Cardolaccia. Data-driven discovery of dimensionless numbers and governing laws: A scaling technique for nuclear test facility design. 2026. cea-05528135v2

HAL Id: cea-05528135

https://cea.hal.science/cea-05528135v2

Preprint submitted on 26 Feb 2026

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 CC BY 4.0 - Attribution - International License
# Data-Driven Discovery of Dimensionless Numbers and Governing Laws: A Scaling Technique for Nuclear Test Facility Design

Camille Razaire,*,a,b Amandine Marrel,a,b Bertrand Iooss,c,d Sébastien Renaudière de Vaux,a and Jérôme Cardolacciaa

<sup>a</sup>CEA, DES, IRESNE
Cadarache, France

<sup>b</sup>Avignon Université, LMA. UPR 2151
Avignon, France

<sup>c</sup>EDF R&D
Chatou, France

<sup>d</sup>SINCLAIR AI Lab
Palaiseau, France

*Email: camille.razaire@cea.fr

Number of pages: 40

Number of tables: 4

Number of figures: 4
Abstract

New designs of water-cooled reactors that include thermal-hydraulics systems undergo safety anal-
ysis during the licensing process. For economic reasons and safety concerns, systems are firstly
tested on reduced scale test facilities. A proper scaling ensures that dominant thermal-hydraulics
safety-related phenomena are captured, even though unavoidable scale distortions occur. Some
existing scaling methods, based on prior knowledge of the physical phenomena at stake, can quan-
tify such distortions. They are however limited when the phenomena are non-linear and coupled,
or even not formalized as an equation. Dimensionless numbers play a central role in scaling
techniques as their scale invariant properties help preserve similarities between reactor and test
facilities. Building on this principle, this paper proposes a data driven alternative to traditional
scaling techniques. From a dataset of physical variables describing the phenomenon of interest, the
method identifies a governing law that captures the dominant safety related phenomenon. This
governing law is expressed in terms of dimensionless numbers, which are physically meaningful
combinations of dimensional variables and are automatically inferred by the algorithm. Based on
a clear mathematical formulation, the proposed data driven method combines advanced regression
analysis with the constrained optimization of a cross validation based objective function. Two vari-
ants are presented: one assuming that the output of the nondimensional governing law is known,
and an extension in which this output is estimated. Both variants identify the dominant input
dimensionless number as well as the explicit form of the governing law. As a proof of concept,
the method is tested on a simulated dataset representative of single-phase natural circulation in a
passive heat removal system.

**Keywords** — Dimensionless Number, Polynomial Regression, Scaling Technique, Passive Safety Systems, Thermal Hydraulics
# I. INTRODUCTION

Improving safety for nuclear reactors is a continuous objective. Hence, new thermal-hydraulics related safety systems are designed and assessed to define their functional limits. However, due to safety and cost requirements, such systems cannot be tested in a full scale nuclear reactor prototype. Therefore, reduced-scale thermal-hydraulics test facilities are used for experiments and licensing process. Appropriate scaling for test facilities design guarantees that safety-relevant phenomena are accurately represented in experiments. Since scaling inherently introduces scale distortions, the scaling analysis requires quantification of these distortions, to scale-up the results from the experiments to the prototypes. From the 1960s, several scaling techniques have been developed for the design of test facilities. Among them, the Power-to-Volume method [1] enabled building full-pressure, full-height and time preserving Integral Test Facilities (ITF) [2], while the Three-Level-Scaling method [3–5] opened the way for ITF with reduced height, area and accelerated time scales [5]. With those techniques, the test facilities designs are based on scaling factors identified and ranked in a PIRT (Phenomena Identification and Ranking Table), built from the analysis of the balance equations. Later, advanced scaling methods emerged with Zuber and colleagues [6, 7], who developed the first technique able to provide higher levels of analysis including distortions quantification: the H2TS method, for Hierarchical-Two-Tiered-Scaling. The implementation of the above-mentioned scaling methods relies on expert knowledge of the governing equations associated with the safety-relevant phenomena. Such dependence was quickly recognized as a major limitation of scaling analyses. Furthermore, advanced scaling methods reach their limits when it comes to quantifying the distortions arising in nonlinear and coupled phenomena.

To mitigate these distortions, scaling factors can be expressed as dimensionless ratios of variables—called dimensionless numbers— that guarantee physical similarity and enhance the understanding of the underlying physical phenomena. Indeed, dimensionless numbers can be interpreted as ratios of characteristic times, forces, or energies, and thus highlight the relative importance and competition between different effects in a same phenomenon [8]. They are widely used in thermal-hydraulics through Vaschy-Buckingham’s theorem [9] that enables their systematic identification. The theorem also shows that a phenomenon can be described with a limited set of nondimensional numbers, rather than the full collection of dimensional variables. As a result, the degrees of freedom of the problem are reduced by the number of fundamental dimensions involved in the initial
physical variables description <sup>a</sup>. However, even after the dimension reduction, the resulting set of nondimensional numbers often remains extensive, and Vaschy-Buckingham's theorem does not provide a ranking of those dimensionless quantities. Thus, the dominant dimensionless numbers linked to the safety-relevant phenomena cannot be identified for scaling purposes.

To tackle the identified limitations of scaling analysis and the aforementioned ranking issue,
data-driven methods are investigated, as they do not require prior knowledge of the phenomena
governing laws and yet enable their prioritization. The data-driven method proposed in this paper
aims to identify the dominant dimensionless number and the associated governing law that will
guide the test facility design, with future work targeting the extension of the approach to two or
more dominant dimensionless numbers. Our method, called D4NL for Data-Driven Discovery of
Dimensionless Numbers and Governing Laws, can thus be considered as a design scaling technique
based on machine learning instead of equation analysis. More generally, it falls within the class of
dimensionless learning methods [10–15], i.e. methods for the learning of dimensionless numbers.

To illustrate this approach, we focus in this paper on the design of passive residual heat re-
moval experiments, whose operating principle relies on single-phase natural circulation. Section II
presents the principles and mathematical formulation of data-driven dimensionless learning with
different state-of-the-art methods. A particular focus is made on Xie et al.’s polynomial regression
method [15] that inspired our D4NL method described in section III. The D4NL method relies on
the minimization of a K-fold-cross-validated loss function using a constrained multistart optimiza-
tion procedure, supported by a space-filling-design strategy tailored to the problem’s constraints.
For this first D4NL variant, the governing law links a known dimensionless output to an unknown
dimensionless input, which is estimated during the learning process. The method is then applied
on a simulated single-phase natural circulation dataset, derived from an analytical model presented
in section IV. The results obtained are presented in section V. In addition, a bootstrap approach is
proposed to estimate confidence intervals on the parameters of the identified dimensionless number.
Finally, Section VI proposes an extension of the D4NL method to the case where the dimensionless
output is unknown. The algorithm estimates it jointly with the input dimensionless number and
the associated governing law. Functional analysis based on variance decomposition and Sobol’
indices are implemented in the algorithm to ease its convergence towards relevant dimensionless

<footnote><sup>a</sup>There are seven fundamental dimensions in physics: L (length), T (time), M (mass), Θ (temperature), I (electric current), N (amount of matter), and J (luminous intensity).</footnote>
output and input.

## II. DATA-DRIVEN LEARNING OF DIMENSIONLESS NUMBERS

### II.A. Principles and mathematical formalization

Dimensionless learning is a data-driven framework designed to uncover governing laws expressed through dimensionless numbers. The central idea is that physical phenomena can often be described more compactly by relating a dimensionless output $\Pi_{\text{out}} \in \mathbb{R}$ to k dimensionless inputs $\{\Pi_{\text{in},(j)}\}_{j=1}^k \in \mathbb{R}^k$. These inputs are constructed from the d physical variables characterizing the phenomenon, $\mathbf{X}_{\text{in}} = (X_1, \dots, X_d)^{\top} \in \mathbb{R}^d$, each carrying physical dimensions. These dimensions can be expressed in terms of $n_{\text{dim,in}}$ fundamental physical dimensions. Formally, the jth dimensionless input is obtained as a monomial combination of the variables in $\mathbf{X}_{\text{in}}$:

$$ \Pi_{\text{in},(j)}(\mathbf{X}_{\text{in}}, \omega_{\text{in}}^{(j)}) = \prod_{i=1}^{d} X_{\text{in},i}^{\omega_{\text{in},i}} \qquad (1) $$

with $\omega_{\text{in}}^{(j)} \in \mathbb{R}^d$ the vector of exponents associated with $\Pi_{\text{in},(j)}$. Each exponent vector $\omega_{\text{in}}^{(j)}$ must satisfy the dimensionless constraint:

$$ \mathbf{A}_{\text{dim,in}} \omega_{\text{in}}^{(j)} = 0 \qquad (2) $$

with $\mathbf{A}_{\text{dim,in}} \in \mathbb{R}^{n_{\text{dim,in}} \times d}$ the fundamental dimension matrix whose columns are the dimension vectors of the physical variables:

$$ \mathbf{A}_{\text{dim,in}} = \begin{bmatrix} \delta(X_1) & \delta(X_2) & \dots & \delta(X_d) \end{bmatrix}. \qquad (3) $$

Here, the mapping function $\delta(\cdot) : \mathbb{R} \to \mathbb{R}^{n_{\text{dim,in}}}$ associates to each physical variable its dimensional vector, i.e. the exponents of the $n_{\text{dim,in}}$ fundamental physical dimensions of the system considered. For instance, in a system where all the variables are expressed in terms of the fundamental dimension basis $(L,T,M)^{\top}$ (i.e. length, time and mass), the dimensional vector of a speed $v$ is $\delta(v) = (1, -1, 0)^{\top}$, corresponding to $L^1T^{-1}M^0$.
The governing relation is then modeled as:

$$ \Pi_{\text{out}} = \mathcal{M}(\Pi_{\text{in}}) + \epsilon_{\text{model}}, \quad (4) $$

where $\mathcal{M}(\cdot): \mathbb{R}^k \to \mathbb{R}$ denotes the predictive model, $\Pi_{\text{in}} = (\Pi_{\text{in},(1)}, \dots, \Pi_{\text{in},(k)})^{\top}$ is the vector with the $k$ dimensionless inputs and $\epsilon_{\text{model}}$ is a random variable representing the discrepancies between the true physical governing law and its approximation by $\mathcal{M}$.

Consider that the sample $\mathcal{D}_N = \{ (\mathbf{X}_{\text{in}}^{(l)}, \Pi_{\text{out}}^{(l)}) \}_{l=1}^N$ is a set of $N$ independent and identically distributed (i.i.d.) measurements, with $\mathbf{X}_{\text{in}}^{(l)} = (X_1^{(l)}, \dots, X_d^{(l)})^{\top} \in \mathbb{R}^d$ the vector of inputs of the $l$th observation, and $\Pi_{\text{out}}^{(l)} \in \mathbb{R}$ the corresponding dimensionless output. Since both inputs and outputs are measured quantities, they are subject to experimental uncertainties. The working relation with data therefore reads

$$ \Pi_{\text{out}}^{(l)} = \mathcal{M} \left( \Pi_{\text{in},(1)}(\mathbf{X}_{\text{in}}^{(l)}, \boldsymbol{\omega}_{\text{in}}^{(1)}), \dots, \Pi_{\text{in},(k)}(\mathbf{X}_{\text{in}}^{(l)}, \boldsymbol{\omega}_{\text{in}}^{(k)}) \right) + \epsilon_{\text{model}} + \epsilon_{\text{meas}}, \quad (5) $$

where $\epsilon_{\text{meas}}$ denotes the measurement noise. For convenience, both error sources are gathered into a single random error term denoted $\epsilon$.

The aim of data-driven dimensionless learning is to infer from $\mathcal{D}_N$ the unknown exponent powers $\{\hat{\boldsymbol{\omega}}_{\text{in}}^{(j)}\}_{j=1}^k$ that define $\{\hat{\Pi}_{\text{in},(j)}(\mathbf{X}_{\text{in}}, \hat{\boldsymbol{\omega}}_{\text{in}}^{(j)})\}_{j=1}^k$, thereby identifying the dominant dimensionless combinations of physical variables. Depending on the methodology used, dimensionless learning can further provide an explicit estimate of the governing law $\mathcal{M}$.

Concerning our context of test facility design, the dimensionless learning method used should respect three requirements:

1. The governing law must be identified in an explicit form and expressed in terms of physically interpretable and meaningful dimensionless numbers.

2. The identification (or estimation) procedure must remain both effective and reliable when applied to scarce datasets, as it is often the case in thermal-hydraulics experiments where the number of available samples is limited.

3. The method must incorporate uncertainty assessment of the estimated exponents, so as to evaluate their statistical identifiability given the available data.
## II.B. State-of-the-art methods

Several state-of-the-art data-driven dimensionless learning methods can be considered to address our problem.

A logarithmic regression is proposed in [10] to estimate $\omega_{in}$ via the following model:

$$ \log(\Pi_{\text{out}}) = \omega_0 + \sum_{i=1}^{d} \omega_{\text{in},i} \log X_{\text{in},i} + \epsilon_L \quad (6) $$

with Gaussian noise $\epsilon_L \sim N(0, \sigma_L^2)$, $\Pi_{out} > 0$ and $X_{in,i} > 0$ for $i \in [1, d]$. This approach restricts the governing law to a single dominant dimensionless number, such as $M = C\Pi_{in,(1)}$ with $C = e^{\omega_0}$. Such a restriction lacks flexibility, as our objective is to generalize the governing law to multiple dominant dimensionless numbers. The logarithmic regression will consequently not be investigated.

Active subspaces used in [11] combine dimension-reduction techniques with global sensitivity analysis to construct unique dimensionless groups. The method is analogous to principal component analysis, except that the projection is performed along directions of largest variation in the model's gradient rather than along dominant variables. Consequently, this method requires substantial amounts of data to reliably estimate the gradient of $M$, and additionally, it does not yield an explicit expression for $M$, thereby failing to meet the prerequisites outlined above.

Neural networks have also been applied to dimensionless learning in [13]. The training procedure minimizes a loss function that includes regularization terms enforcing both sparsity of $\omega_{in}$ and compliance with the dimensionless constraint (cf. Eq. (2)). As with active subspaces, this method does not provide an explicit expression of the model $M$. Moreover, the simultaneous estimation of $\omega_{in}$ and $M$ proves challenging in the small-data context typical of experimental campaigns on nuclear system test facilities.

The polynomial regression algorithm introduced by Xie et al. in [15] estimates the governing law by minimizing a loss function between $\Pi_{out}$ and a polynomial model of degree $P$ built on $k$ dominant dimensionless inputs $\Pi_{in}$, each one parametrized by a couple $(\gamma_{in}, \beta)$. The parameters $\gamma_{in} \in \mathbb{R}^{d-n_{\text{dim,in}}}$ are linked to $\omega_{in} \in \mathbb{R}^d$ through the dimensionless constraint (more details in Section II.C), thereby fixing the model form, while the regression coefficients $\beta \in \mathbb{R}^{P+1}$ adjust the model to the data. In the examples treated in [15], $P \le 5$ and the polynomial model involves
one or two dominant dimensionless numbers. The optimal parameters ($\hat{\gamma}_{in}, \hat{\beta}$) are obtained via a two-level optimization scheme, with one loop for each group of parameters: a grid search for the structural exponents $\gamma_{in}$, and a gradient-based optimization for the regression coefficients $\beta$. This hierarchical separation acts as an implicit regularization, limiting compensatory effects between the model's structure and its fit to the data. It also reduces overfitting risks in small-data regimes. This method provides an explicit form of the estimated governing law and can be applied to small datasets without excessive risk of overfitting, provided the polynomial degree remains moderate. Furthermore, the polynomial regression model offers greater flexibility than logarithmic regression, which is restricted to a single dominant dimensionless number. By allowing multiple dimensionless inputs and nonlinear interactions through polynomial terms, the approach can capture a wider range of governing laws. However, this increased expressiveness must be balanced against the risk of overfitting, which can be mitigated through additional regularization techniques. Xie et al.'s method will therefore serve as a starting point for the development of D4NL method.

Nevertheless, some drawbacks from Xie et al.'s method [15] have been identified. These include the suboptimal two-level iterative optimization scheme and the lack of uncertainty assessment, either with respect to model-form assumptions or to measurement data. Furthermore, no uncertainty quantification is performed in the analysis of the results. To facilitate generalization, we therefore propose improvements such as a more rigorous mathematical formalization of the model hypotheses and the adoption of advanced statistical tools (e.g. K-Fold cross-validation, a posteriori uncertainty inference). These enhancements pave the way for extended dimensionless learning, that could identify nonlinear governing laws, with multiple dimensionless numbers, or even with unknown dimensionless output $\Pi_{out}$. Before presenting in detail our improved version of Xie et al.'s method, named D4NL, we briefly recall the problem reparametrization introduced in [15] that we also adopt.

## II.C. Preliminary problem reparametrization

As explained in [15], searching for $\omega_{in}^{(j)}$ that satisfies Eq. (2) is equivalent to restricting $\omega_{in}^{(j)}$ to the null space of $A_{dim,in}$. In other words, $\omega_{in}^{(j)}$ can be expressed as:

$$ \omega_{in}^{(j)} = V_{in} \gamma_{in}^{(j)} \qquad (7) $$
with:

- $V_{in} \in \mathbb{R}^{d \times (d-n_{dim,in})}$ an orthonormal basis of the null space of $A_{dim,in}$, obtained via the Gram-Schmidt process [16] and referred to as the dimensionless basis.

- $\gamma_{in}^{(j)} \in D_{\gamma_{in}} \subset \mathbb{R}^{(d-n_{dim,in})}$ the dimensionless basis coefficients of the $j$th dimensionless number.

This reparametrization reduces the problem from $d$ degrees of freedom, based on the choice of $\omega_{in}^{(j)}$, to $d - n_{dim,in}$ degrees of freedom, where the task is to identify the coefficients $\gamma_{in}^{(j)}$. The column ordering in $A_{dim,in}$ is arbitrary. Therefore, the basis vectors composing $V_{in}$ are not unique. However, any valid choice of $V_{in}$ spans the same null space and leads to equivalent results.

As the problem is now $\gamma_{in}$-dependent, the definition of $\Pi_{in,(j)}$ in Eq. (1) becomes

$$ \Pi_{in,(j)}(\mathbf{X}_{in}, \omega_{in}^{(j)}) = \prod_{i=1}^{d} X_{in,i}^{\omega_{in,i}^{(j)}} = \prod_{i=1}^{d} X_{in,i}^{(V\gamma)_{in,i}^{(j)}} = \Pi_{in,(j)}(\mathbf{X}_{in}, \gamma_{in}^{(j)}) \quad (8) $$

and the dimensionless learning method now aims at finding the optimal $\hat{\gamma}_{in}^{(j)}$ that generates $\hat{\Pi}_{in,(j)}(\mathbf{X}_{in}, \hat{\gamma}_{in}^{(j)})$.

### III. D4NL METHOD: DATA-DRIVEN DISCOVERY OF DIMENSIONLESS NUMBERS AND GOVERNING LAWS

As an initial step toward dimensionless learning, the development of the model relies on the following simplifying assumptions—which will be progressively relaxed in future extensions:

- The governing law depends on a single dominant dimensionless number ($k=1$). Hence, exponent ($j$) is omitted and $(\omega_{in}^{(j)}, \gamma_{in}^{(j)})$ becomes $(\omega_{in}, \gamma_{in})$.

- The governing law is a polynomial model function of this input dimensionless number and is linear in its parameters $\beta$. Identifying the law reduces to a linear estimation problem.

- The dimensionless output $\Pi_{out}$ is assumed to be known.

Under these assumptions, the following sections present the construction of the Data-Driven Discovery of Dimensionless Numbers and governing Laws (D4NL) method.
### III.A. Polynomial regression for one dominant input dimensionless number

Many physical correlations between dimensionless numbers are expressed as power laws or polynomial relations. Assuming a polynomial form facilitates the identification of the model parameters and provides its explicit analytical form. Following the approach of Xie et al. [15], we propose a dimensionless learning process based on polynomial regression, such that:

$$ \Pi_{\text{out}} = f(\beta, \Pi_{\text{in}}(\mathbf{X}_{\text{in}}, \gamma_{\text{in}})) + \epsilon \quad (9) $$

with

$$ f(\beta, \Pi_{\text{in}}(\mathbf{X}_{\text{in}}, \gamma_{\text{in}})) = f(\mathbf{X}_{\text{in}}, \beta, \gamma_{\text{in}}) = \sum_{l=0}^{P} \beta_l (\Pi_{\text{in}}(\mathbf{X}_{\text{in}}, \gamma_{\text{in}}))^l \quad (10) $$

where $f$ denotes the polynomial model, $\beta = (\beta_0, \beta_1, \dots, \beta_P)^\top \in \mathbb{R}^{P+1}$ are the regression coefficients, and $P$ is kept limited to mitigate overfitting (see Sections III.C and V). Alternatively, penalized regression methods such as Ridge or LASSO can also be employed [17].

The polynomial model assumption may induce some modeling error if the true physical relationship deviates from a polynomial form. In such cases, $\mathcal{M}$ captures only part of the phenomenon, while the unmodeled effects are included into the error term $\epsilon$, together with measurement uncertainties. Regarding the statistical nature of $\epsilon$, note that measurement errors on physical quantities propagate multiplicatively when forming dimensionless numbers, which would typically lead to log-normal noise distributions. However, in this work we assume a global Gaussian error model, $\epsilon \sim \mathcal{N}(0, \sigma_{\epsilon}^2)$. This assumption is motivated by the fact that the dominant source of error is expected to be the modeling approximation (due to the simplicity of the polynomial form), rather than the measurement noise. Under this reasoning, the Gaussian hypothesis provides a rational and tractable approximation. This assumption should nevertheless be validated in practice through residual analysis (see Section V.B).

In the next sections, we turn to the practical aspects of the D4NL method. Having established the polynomial form of the governing law, the learning task therefore reduces to the estimation of the parameters $\beta$ and $\gamma_{\text{in}}$ from the dataset $\mathcal{D}_N$. This amounts to identifying their optimal values $\hat{\beta}$ and $\hat{\gamma}_{\text{in}}$, which define the learned model $\hat{f}(\mathbf{X}_{\text{in}}) := f(\mathbf{X}_{\text{in}}, \hat{\beta}, \hat{\gamma}_{\text{in}})$. We therefore introduce in the following sections the loss function, as well as the optimization algorithm employed to minimize it and to estimate the parameters from $\mathcal{D}_N$.
Hereafter, the observed dataset is denoted in the following $\{(\mathbf{x}_{\text{in}}^{(j)}, \pi_{\text{out}}^{(j)})\}_{j=1}^N$, where low- ercase notations refer to observed realizations, while uppercase notations are reserved for random variables. This convention does not apply to the physical variables introduced in Section IV.

### III.B. Cross-validation-based loss function for parameter estimation

The construction of the loss function in D4NL relies on two key principles. First, to ensure robustness against small datasets and to mitigate overfitting, the measure of fit to the data is computed through K-Fold Cross-Validation (CV). By repeatedly splitting the dataset into training and validation folds, CV provides an unbiased estimate of the generalization error while limiting the loss of information that would occur with a single train/test split [17, 18]. Second, for a fixed choice of $\gamma_{\text{in}}$ (i.e. for a given model structure), the regression coefficients $\beta$ admit a closed-form least squares solution. This property allows the estimation of $\beta$ to be embedded within the CV procedure, so that the loss function depends only on $\gamma_{\text{in}}$ and can be minimized accordingly.

Formally, from a given n-size learning sample $\mathcal{D}_n = \{(\mathbf{x}_{\text{in}}^{(j)}, \pi_{\text{out}}^{(j)})\}_{j=1}^n$, we define the vector of observed outputs $\pi_{\text{out}} := (\pi_{\text{out}}^{(1)}, \pi_{\text{out}}^{(2)}, \dots, \pi_{\text{out}}^{(n)})^{\top}$ and the vector of predictions $\hat{f}_{\text{CV}}(\gamma_{\text{in}})$ obtained by K-Fold CV for a given $\gamma_{\text{in}}$:

$$ \hat{f}_{\text{CV}}(\gamma_{\text{in}}) := \left\{ \hat{f} \left( \mathbf{x}_{\text{in}}^{(\{k\})}, \hat{\beta}_{-\{k\}}(\gamma_{\text{in}}), \gamma_{\text{in}} \right) \right\}_{k=1}^K \quad (11) $$

where $\{k\}$ denotes the validation sample of the k-th fold, and $\hat{\beta}_{-\{k\}}(\gamma_{\text{in}})$ is the least squares estimator of $\beta$ computed on all data except the validation fold. The optimal set of parameters $\hat{\gamma}_{\text{in}}(\mathcal{D}_n)$ is then obtained as:

$$ \hat{\gamma}_{\text{in}}(\mathcal{D}_n) = \underset{\gamma_{\text{in}} \in \mathcal{D}_{\gamma_{\text{in}}}}{\operatorname{argmin}} \left( \mathcal{L}\left(\pi_{\text{out}}, \hat{f}_{\text{CV}}(\gamma_{\text{in}})\right) \right), \quad (12) $$

where $\mathcal{D}_{\gamma_{\text{in}}}$ is the feasible domain of $\gamma_{\text{in}}$ (defined in the next section) and $\mathcal{L}$ denotes the loss function. $\mathcal{L}$ is chosen here as the Mean Square Error (MSE)

$$ \mathcal{L}\left(\pi_{\text{out}}, \hat{f}_{\text{CV}}(\gamma_{\text{in}})\right) = \frac{1}{n} \sum_{j=1}^{n} \left( \pi_{\text{out}}^{(j)} - \hat{f}_{\text{CV}}(\mathbf{x}_{\text{in}}^{(j)}, \gamma_{\text{in}}) \right)^2 . \quad (13) $$

Thanks to the K-Fold CV procedure, the model predictions are always compared to the observed
values $\pi_{\text{out}}$ on data unseen during the estimation of $\beta$, ensuring that the loss function provides an unbiased measure of the generalization performance. In this nested scheme, $\beta$ is treated as an inner parameter estimated analytically for each candidate $\gamma_{\text{in}}$, while $\gamma_{\text{in}}$ is optimized externally via cross-validation.

Once $\hat{\gamma}_{\text{in}}(\mathcal{D}_n)$ is identified via (12) (and the solving of the optimization problem explained in Section III.C), the corresponding coefficients $\hat{\beta}(\hat{\gamma}_{\text{in}}(\mathcal{D}_n))$ are obtained by fitting the linear regression on the entire dataset $\mathcal{D}_n$. The final governing law estimator can then be applied to any new input $\mathbf{x}_{\text{in}}^*$, and is written as:

$$ \hat{f}(\mathbf{x}_{\text{in}}^*, \hat{\gamma}_{\text{in}}(\mathcal{D}_n)) := f(\mathbf{x}_{\text{in}}^*, \hat{\beta}(\hat{\gamma}_{\text{in}}(\mathcal{D}_n)), \hat{\gamma}_{\text{in}}(\mathcal{D}_n)), \quad (14) $$

which defines the model predictions for unseen data as:

$$ \hat{\pi}_{\text{out}}(\mathbf{x}_{\text{in}}^*, \hat{\gamma}_{\text{in}}(\mathcal{D}_n)) = \hat{f}(\mathbf{x}_{\text{in}}^*, \hat{\gamma}_{\text{in}}(\mathcal{D}_n)). \quad (15) $$

### III.C. Multi-start optimization algorithm in convex polytope parameter space

As introduced in Eq. (12), the optimal set of parameters $\hat{\gamma}_{\text{in}}(\mathcal{D}_n)$ is obtained by minimizing the loss function over the admissible domain $\mathcal{D}_{\gamma_{\text{in}}}$. This domain is derived from the physical constraints imposed on $\omega_{\text{in}}$, in addition to the nondimensional requirement. To guarantee physically meaningful solutions, each coordinate $\omega_{\text{in},i}$ ($i \in [1, d]$) is restricted to lie within $[-\omega_{\text{in}}^{\text{max}}, \omega_{\text{in}}^{\text{max}}]$. Although sparsity of $\omega_{\text{in}}$ could also be enforced through regularization, this aspect will not be considered here. A trade-off must be found between the choice of $\omega_{\text{in}}^{\text{max}}$ and the polynomial degree $P$: too large values of $\omega_{\text{in},i}$ ($i \in [1, d]$) reduce the physical interpretability of the generated $\pi_{\text{in}}(\mathbf{x}_{\text{in}}, \omega_{\text{in}})$, while too small values increase the required polynomial degree $P$ of the constructed governing law, which may lead to overfitting. In the following, we set $\omega_{\text{in}}^{\text{max}} = 3$ as a reasonable and general constraint, since most dimensionless numbers reported in the literature do not contain variables raised to powers higher than 3.

The constraint applied on the $\omega_{\text{in},i}$ ($i \in [1, d]$) values directly impacts the feasible domain $\mathcal{D}_{\gamma_{\text{in}}}$ of $\gamma_{\text{in}}$. Since $\gamma_{\text{in}} = \mathbf{V}_{\text{in}}^{\top}\omega_{\text{in}}$, the feasible set is the image of the hypercube $[-\omega_{\text{in}}^{\text{max}}, \omega_{\text{in}}^{\text{max}}]^d$ under the linear transformation $\mathbf{V}_{\text{in}}^{\top}$. $\mathcal{D}_{\gamma_{\text{in}}}$ is therefore a convex polytope, delimited by hyperplanes
inherited from the bounds on $\omega_{in}$:

$$ \mathcal{D}_{\gamma_{in}} := \{ \gamma_{in} \mid \mathbf{V}_{in} \gamma_{in} \in [-\omega_{in}^{\max}, \omega_{in}^{\max}]^d \}. \quad (16) $$

An illustration of such feasible domain for a problem involving $d = 8$ variables described by $n_{dim,in} = 4$ fundamental physical dimension (yielding in $\gamma_{in} \in \mathbb{R}^4$) is provided in Figure 1. The blue points in Figure 1 correspond to admissible realizations of $\gamma_{in}$ for a given matrix $\mathbf{V}_{in}$ and $\omega_{in}^{\max} = 3$.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Example scatter plots of admissible realizations of $\gamma_{in}$ (in blue) and selected support points (in orange), projected onto each pairs of coordinate.</figcaption></figure>

$\mathcal{D}_{\gamma_{in}}$ defines the admissible search space in which the optimization of the loss function (Eq. (12)) must be carried out. To solve this constrained minimization problem, a gradient-descent based algorithm (e.g., quasi-Newton type such as Sequential Least Squares Programming (SLSQP) Algorithm) is employed, specifically designed to handle inequality constraints defining the feasible domain $\mathcal{D}_{\gamma_{in}}$. The approach is reinforced by a multi-start strategy: several indepen-
dent optimizations are launched from distinct initial points to explore $D_{\gamma_{in}}$ thoroughly and to mitigate the risk of the descent procedure being trapped in local minima. In order to maximize the representativity of these starting points, they must uniformly cover $D_{\gamma_{in}}$. For this purpose, advanced space-filling techniques are employed, and more specifically the iterative support points method [19]. This method constructs a distributionally representative subset of points (orange dots in Figure 1, noted $N_{SP,in}$) from a large pool of candidates sampled within $D_{\gamma_{in}}$. Compared to a purely random selection, support points provide a more uniform and informative coverage of the convex polytope, thereby improving the likelihood of reaching the global minimizer.

This combination of constrained gradient descent with multi-start initialization provides a pragmatic balance between robustness against local minima and computational efficiency compared to more complex global optimization methods. It is computationally efficient enough to be embedded within a nested cross-validation framework for the assessment of the final model. The algorithm of the D4NL estimation is summarized by Algorithm 1 in Appendix A.

### III.D. Statistical metrics and diagnostics for final validation

Once the optimal parameters $\hat{\gamma}_{in}(D_n)$ have been estimated, the predictive quality of the resulting estimator $\hat{\pi}_{out}(x_{in}, \hat{\gamma}_{in}(D_n))$ must be assessed on an unseen dataset (test sample of size $n_T$). For this purpose, we rely on the predictivity coefficient $Q^2$, a standard metric in statistical modeling that quantifies the proportion of variance explained in prediction. Formally, for a test dataset of size $n_T$, $Q^2$ is defined as:

$$Q^2 = 1 - \frac{\sum_{j=1}^{n_T} \left( \pi_{out}^{(j)} - \hat{\pi}_{out}^{(j)} \right)^2}{\sum_{j=1}^{n_T} \left( \pi_{out}^{(j)} - \bar{\pi}_{out} \right)^2} = 1 - \frac{\text{PRESS}}{\text{TSS}}, \quad (17)$$

where $\pi_{out}^{(j)}$ and $\hat{\pi}_{out}^{(j)}$ respectively denote the observed output and the model prediction for the jth observation, and $\bar{\pi}_{out}$ is the mean of the observed outputs in the test dataset. PRESS is the Predictive Residual Error Sum of Squares, and TSS the Total Sum of Squares [17]. The interpretation of $Q^2$ is straightforward: the closer $Q^2$ is to 1, the better the predictive accuracy (see [20] for more details on $Q^2$). This metric can also be supplemented by graphical analyses of the residuals $\{\pi_{out}^{(j)} - \hat{\pi}_{out}^{(j)}\}_{j=1}^{n_T}$.

In practice, residuals and the associated $Q^2$ are computed on a test dataset. However,
when the total sample size is limited, a simple train/test split is not feasible. Following the same
rationale as in Section III.B, these quantities can instead be estimated using an outer K'-Fold
Cross-Validation scheme. This outer loop provides an unbiased evaluation of the final model, while
the inner loop (based on K-Fold CV) is already used for the optimization of $\gamma_{in}$ (as described in
Eq. (12)). Formally, for each outer validation fold $k' \in \{1, \dots, K'\}$, the optimization algorithm is
run on the reduced learning sample $\mathcal{D}_n \setminus \{k'\}$, excluding the data of the $k'$-th outer fold. Within this
reduced dataset, the inner K-Fold CV is applied to estimate $\hat{\gamma}_{in}(\mathcal{D}_{n \setminus \{k'\}})$. The resulting estimator
is then used to predict the outputs of the left-out fold $k'$, yielding residuals $\pi_{\text{out}}^{(j)} - \hat{\pi}_{\text{out}}^{(j)}$ for all $j \in k'$.
Aggregating predictions across all $K'$ folds produces the cross-validated prediction vector $\hat{f}_{\text{CV,outer}}$,
from which the final $Q_{\text{CV,outer}}^2$ is computed. This nested procedure ensures that the reported $Q^2$
reflects the true generalization ability of the governing law estimator, since model selection (inner
loop) and model validation (outer loop) are strictly separated. The entire procedure is summarized
by Algorithm 2 in Appendix A.

The following paragraph presents the first application of the D4NL method to a single-phase
natural circulation dataset for the design of a passive heat removal system.

## IV. DESCRIPTION OF THE SINGLE-PHASE THERMOSYPHON ANALYTICAL TEST CASE

Thermosyphons, or natural circulation systems, are passive heat transport devices. In the
nuclear field, they are proposed for passive core cooling systems applications and tested in sev-
eral experimental loops around the world (e.g., PASI [21], CIRCUS [22], EXOCET [23]). Ther-
mosyphons are also an active research topic, and different mechanisms of instabilities are described
in the literature, both for single-phase [24] and two-phase flows [21–23, 25–28]. The main results
and observations of natural circulation are summarised in Vijayan’s book [29]. In particular, for
single-phase natural circulation, the nondimensional fluid velocity (i.e. *Re* number) depends on
few known dimensionless numbers (see Eq. (26)). Hence, single-phase natural circulation is a good
candidate for validating the D4NL method. The test facility studied is a closed single-phase ther-
mosyphon with uniform diameter and adiabatic parts between the horizontals heat source and sink
(see Figure 2 for the system’s sketch). Section IV.A presents the steady state model of single-phase
natural circulation established in the loop. The model is then used in section IV.B to generate the
data needed by the dimensionless learning method.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Sketch of a thermosyphon with horizontals heat source (red rectangular) and heat sink (blue rectangular).</figcaption></figure>

IV.A. Physical model

IV.A.1. *Dimensional formulation*

Natural circulation flow results from the balance between buoyancy and pressure losses. The flow is driven by the gravitational head, created by a vertical thermally-induced density gradient between a heat sink located at an elevation *H* above a heat source. To obtain the steady state behavior of such flow, several approximations are made:

- Flow evolution is considered one-dimensional. Thereby, the natural circulation bulk velocity is averaged over the whole loop.

- The analysis is limited to single-phase liquid conditions for water. We thus respect newtonian and incompressible conditions.

- Boussinesq's approximation is considered to be valid [30]. Consequently, the fluid properties in the loop are assumed constant, except for the density in the buoyancy force term which
varies linearly with the temperature:

$$ \rho(T) = \rho_0 [1 - \beta_{th}(T - T_0)] \quad (18) $$

with T the temperature, $\rho_0$ the homogeneous fluid density at average temperature $T_0$ and $\beta_{th}$ the thermal expansion coefficient at $T_0$.

- Local losses are neglected compared to frictional losses. Physically, this would correspond to a loop having long pipes with smooth elbows.

Under the previous assumptions, integrating the momentum conservation equations over the entire loop yields the following 0-D model:

$$ \Delta P_g + \Delta P_f = 0 \quad (19) $$

with $\Delta P_g$ the gravitational head and $\Delta P_f$ the frictional pressure loss, equal to:

$$ \Delta P_g = \Delta \rho g H = -\rho_0 g H \beta_{th} \Delta T \quad (20) $$

$$ \Delta P_f = f_r \frac{\rho_0 v^2}{2} \frac{L_t}{D} \quad (21) $$

with $f_r$ the friction coefficient, $v$ the average fluid velocity, $L_t$ the total piping length (horizontal and vertical sections), $D$ the pipe diameter, $g$ the gravity acceleration and $\Delta T = T_h - T_c > 0$ the temperature difference between the hot and the cold legs in the loop.

By substituting Eqs. (21) and (20) in Eq. (19), we readily obtain the steady-state flow velocity $v$ for single-phase natural circulation

$$ v = \sqrt{\frac{2HD\beta_{th}\Delta Tg}{f_rL_t}} \quad (22) $$

### IV.A.2. Nondimensional formulation

The previous fluid velocity equation is nondimensionalized using the parameters ratio $\frac{\mu}{\rho_0 D}$ (homogeneous to a velocity) with $\mu$ the average dynamic viscosity, such that:

$$ \frac{\rho_0 v D}{\mu} = \sqrt{\frac{2HD^3 \rho_0^2 \beta_{th} \Delta T g}{f_r L_t \mu^2}} \quad (23) $$
We therefore identify two dimensionless numbers in Eq. (23):

- The Reynolds number (Re) defined as:

$$Re = \frac{\rho_0 v D}{\mu} \qquad (24)$$

that quantifies the competition in a flow between inertial forces (at the numerator) and viscous forces (at the denominator).

- The modified Grashof number (Gr_m) given by:

$$Gr_m = \frac{D^3 \rho_0^2 \beta_{th} \Delta T g}{\mu^2} \qquad (25)$$

which characterizes natural convection as the ratio of the “characteristic times for the diffusion of momentum and the free fall time” (cf. Lykoudis in [8]).

Hence, Eq. (23) can be expressed as follows:

$$Re = \left( \frac{2H}{f_r L_t} Gr_m \right)^{1/2} \qquad (26)$$

In circular pipes, the friction coefficient $f_r$ is analytically defined for laminar regime, whereas Blasius' empirical correlation [31] is often used for turbulent regime. In this approach, the laminar-turbulent transition is not modeled, while laminar and turbulent expressions are defined as:

$$f_r = \frac{a}{Re^b} \quad \text{with} \quad \begin{cases} a = 64, b = 1 & \text{if } Re < 2000 \\ a = 0.3164, b = 0.25 & \text{otherwise.} \end{cases} \qquad (27)$$

With Eq. (27), Eq. (26) becomes

$$Re = \begin{cases} \frac{1}{32} \frac{H}{L_t} Gr_m & \text{for } Re \le 2000, \\ \left( \frac{2H}{0.3164 L_t} Gr_m \right)^{4/7} & \text{otherwise.} \end{cases} \qquad (28)$$
By approximating the numerical constants and introducing $\tilde{Gr}_m = \frac{H}{L_t}Gr_m$, we obtain:

$$ Re = \begin{cases} 0.0313 \tilde{Gr}_m & \text{for } Re \le 2000 \\ 2.87 (\tilde{Gr}_m)^{4/7} & \text{otherwise} \end{cases} \qquad (29) $$

which defines two governing laws linking two dimensionless numbers, $Re$ and $\tilde{Gr}_m$.

In the following, the aim will be to automatically recover the expression of $\tilde{Gr}_m$ as a function of dimensional inputs, as well as the governing law that links it to the dimensionless output $Re$.

## IV.B. Simulated dataset

For the single-phase thermosyphon test case, the dimensional input vector is defined as $\mathbf{X} = (D, \mu, \rho_0, \beta_{th}, \Delta T, H, L_t, g)^{\top}$ (then $d = 8$) and the dimensionless output is the Reynolds number:

$$ \Pi_{\text{out}} := Re = f(D, \mu, \rho_0, \beta_{\text{th}}, \Delta T, H, L_t, g). \qquad (30) $$

Although $g$ does not vary in practice, it is treated as an input variable in the dimensionless learning process. The problem is split between laminar ($Re \le 2000$) and turbulent regimes ($Re > 2000$): distinct expressions for the true governing law $f$ are given by Eq. (29) for each regime. Each input variable $X_i$ ($i \in [1, d]$) is sampled from a uniform distribution $U([a_i, b_i])$, whose bounds $[a_i, b_i]$ are specific to each regime and specified in Table I.

Based on these definitions, two synthetic datasets $\{ (\mathbf{x}_{\text{in}}^{(j)}, \pi_{\text{out}}^{(j)}) \}_{j=1}^N$ are generated, one for each regime, by sampling $N = 200$ independent and identically distributed realizations of the input vector $\mathbf{x}_{\text{in}}$ and computing the associated outputs $\pi_{\text{out}}$ with Eq. (29). This choice of sample size is consistent with the order of magnitude of datasets obtained in experimental campaigns on thermosyphons. Moreover, to reproduce realistic conditions, Gaussian noise $\epsilon \sim N(0, \sigma_{\epsilon})$ is added to the outputs. This perturbation accounts both for model error—i.e. the discrepancy between the true governing law and the polynomial form assumed in the D4NL method—and for measurement uncertainties on $\pi_{\text{out}}$. The noise level is scaled to the mean Reynolds number of the dataset ($\sigma_{\epsilon} = \tau \overline{Re}$, with $\tau \in \mathbb{R}^{+*}$). The resulting datasets are denoted generically as $D_{200,\text{lam}}(\tau)$ and $D_{200,\text{tur}}(\tau)$, where each dataset corresponds to $N = 200$ samples generated for a given noise level $\tau$.
<table><caption>TABLE I
Physical inputs of the thermosyphon model and associated probability distributions.</caption><thead><tr><th rowspan="2">Input variable</th><th rowspan="2"></th><th rowspan="2">Notation</th><th colspan="2">[a, b] boundaries of U(a, b)</th><th rowspan="2">Unit</th><th rowspan="2">Dimension</th></tr><tr><th>Laminar regime</th><th>Turbulent regime</th></tr></thead><tbody><tr><td>X<sub>1</sub></td><td>Pipe diameter</td><td>D</td><td>[9; 10] × 10<sup>-3</sup></td><td>[40; 55] × 10<sup>-3</sup></td><td>m</td><td>L<sup>1</sup></td></tr><tr><td>X<sub>2</sub></td><td>Viscosity</td><td>μ</td><td>[6.2; 7.5] × 10<sup>-4</sup></td><td>[5.0; 7.5] × 10<sup>-4</sup></td><td>Pa.s</td><td>L<sup>-1</sup>T<sup>-1</sup>M<sup>1</sup></td></tr><tr><td>X<sub>3</sub></td><td>Density</td><td>ρ<sub>0</sub></td><td>[991; 995]</td><td>[985; 994]</td><td>kg.m<sup>-3</sup></td><td>L<sup>-3</sup>M<sup>1</sup></td></tr><tr><td>X<sub>4</sub></td><td>Thermal expansion coefficient</td><td>β<sub>th</sub></td><td>[3.2; 4.0] × 10<sup>-4</sup></td><td>[3.2; 4.5] × 10<sup>-4</sup></td><td>K<sup>-1</sup></td><td>Θ<sup>-1</sup></td></tr><tr><td>X<sub>5</sub></td><td>Temperature difference</td><td>ΔT</td><td>[5; 20]</td><td>[10; 50]</td><td>°C</td><td>Θ<sup>1</sup></td></tr><tr><td>X<sub>6</sub></td><td>Loop height</td><td>H</td><td>[1.5; 3]</td><td>[1.5; 3]</td><td>m</td><td>L<sup>1</sup></td></tr><tr><td>X<sub>7</sub></td><td>Total loop length</td><td>L<sub>t</sub></td><td>[3; 18]</td><td>[3; 18]</td><td>m</td><td>L<sup>1</sup></td></tr><tr><td colspan="2">Constant input</td><td>Notation</td><td colspan="2">Value</td><td>Unit</td><td>Dimension</td></tr><tr><td>X<sub>8</sub></td><td>Gravity acceleration</td><td>g</td><td colspan="2">9.81</td><td>m.s<sup>-2</sup></td><td>L<sup>1</sup>T<sup>-2</sup></td></tr></tbody></table>

In this test case, from the considered input physical variables, the fundamental dimension matrix associated to Π<sub>in</sub> for nondimensional constraint is:

$$ \mathbf{A}_{\text{dim,in}} = \begin{pmatrix} D & \mu & \rho_0 & \beta_{th} & \Delta T & H & L_t & g \\ 1 & -1 & -3 & 0 & 0 & 1 & 1 & 1 \\ 0 & -1 & 0 & 0 & 0 & 0 & 0 & -2 \\ 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 & 1 & 0 & 0 & 0 \end{pmatrix} \begin{array}{l} L \\ T \\ M \\ \Theta \end{array} \quad (31) $$

## V. APPLICATION OF THE D4NL METHOD TO THE THERMOSYPHON ANALYTICAL TEST CASE WITH Π<sub>out</sub> KNOWN

The D4NL method is applied to the thermosyphon datasets in order to assess its ability to rediscover the known governing laws under controlled yet imperfect conditions. In this setting, the governing laws to be identified is known a priori and corresponds to a first-degree polynomial (P = 1 and β<sub>0</sub> = 0 in Eq. (10)). This is a favorable configuration, fully consistent with the polynomial hypothesis assumed for the functional form of f in Eq. (10). In practice, for unknown governing laws, it is advised to fix the polynomial degree to a low value, and progressively increase
it if residuals exhibit non-Gaussian behavior (see Section V.B on residuals analysis) or if the identified dimensionless number exponents are too high.

Datasets are generated with a noise variance parameter $\tau$ varied between 0.05 and 0.2, enabling the evaluation of the method's robustness to increasing levels of noise contamination. The D4NL algorithm (Algo 2 in Appendix A) is applied assuming $P = 1$ and with $K = K' = 8$ folds in the nested cross-validation procedure and $N_{SP,in} = 100$ for the multistart process. This proof-of-concept experiment therefore provides a stringent test of the algorithm's ability to correctly identify the dominant dimensionless number and recover the governing law under favorable but noisy conditions. In subsequent work, the assumption of $P = 1$ could be relaxed to polynomials of higher degree, allowing the evaluation of the method's capacity to capture more complex relationships and to infer the governing law in its entirety.

## V.A. Results

The results obtained with the D4NL method, applied to $N = 200$-sample datasets contaminated by different noise levels for both regimes, are summarized in Table II and illustrated in Figure 3. The exponents inferred by the algorithm in Table II are reported as rational numbers, within a tolerance of 0.1. Results illustrate the method's ability to correctly identify $\tilde{Gr}_m$ and relate it to the Reynolds number via a monomial form. For the laminar regime, the recovered expressions are in perfect agreement with the theoretical solution. In the turbulent regime, slight deviations appear at higher noise levels. Remarkably, even when the noise level is relatively high (standard deviation up to 20% of the mean Reynolds number), the correct dimensionless formulation is consistently recovered. The noise impact is observed primarily on the governing law parameter $\beta$, while the structure of the dimensionless number remains stable. In terms of predictive performance, the $Q_{CV}^2$ values reported in Table II equals $Q_{max}^2$, the theoretical maximum achievable under noisy conditions. This upper bound corresponds to 1 minus the fraction of output variance attributable to the added noise, i.e., the part of the signal that no model can recover. Figure 3 further confirms this agreement: simulated noisy values of $Re$ plotted against the estimated dimensionless numbers align closely with the linear model $\hat{\beta}(\hat{\gamma}_{in})\hat{\pi}_{in}(\mathbf{x}_{in}, \hat{\gamma}_{in})$ (solid red line).

It was observed that less starting points in the multistart process ($N_{SP,in} \approx 20$) also yields satisfactory results, without loss of information.
<caption>TABLE II</caption>

Estimated values of $\hat{\beta}(\hat{\gamma}_{in})$, $\hat{\pi}_{in}(\mathbf{x}, \hat{\gamma}_{in})$ and the predictive performance metrics for different noise levels.

<table><thead><tr><th rowspan="2">Noise standard<br>deviation level<br>τ (σ<sub>e</sub> = τR<sub>e</sub>)</th><th colspan="4">Laminar regime</th><th colspan="4">Turbulent regime</th></tr><tr><th>β̂(γ̂<sub>in</sub>)</th><th>π̂<sub>in</sub>(x<sub>in</sub>, γ̂<sub>in</sub>)</th><th>Q<sup>2</sup><sub>CV</sub></th><th>Q<sup>2</sup><sub>max</sub></th><th>β̂(γ̂<sub>in</sub>)</th><th>π̂<sub>in</sub>(x<sub>in</sub>, γ̂<sub>in</sub>)</th><th>Q<sup>2</sup><sub>CV</sub></th><th>Q<sup>2</sup><sub>max</sub></th></tr></thead><tbody><tr><td>0.05</td><td>0.029</td><td>D<sup>3</sup>ρ<sub>0</sub><sup>2</sup>β<sub>th</sub>ΔTHg / μ<sup>2</sup>L<sub>t</sub></td><td>0.99</td><td>0.99</td><td>2.7</td><td>(D<sup>3</sup>ρ<sub>0</sub><sup>2</sup>β<sub>th</sub>ΔTHg / μ<sup>2</sup>L<sub>t</sub>)<sup>4/7</sup></td><td>0.98</td><td>0.98</td></tr><tr><td>0.1</td><td>0.027</td><td>D<sup>3</sup>ρ<sub>0</sub><sup>2</sup>β<sub>th</sub>ΔTHg / μ<sup>2</sup>L<sub>t</sub></td><td>0.97</td><td>0.97</td><td>2.5</td><td>(D<sup>29</sup>ρ<sub>10</sub><sup>2</sup>β<sub>th</sub>ΔTH<sup>11</sup>g / μ<sup>2</sup>L<sub>t</sub>)<sup>4/7</sup></td><td>0.94</td><td>0.94</td></tr><tr><td>0.2</td><td>0.023</td><td>D<sup>3</sup>ρ<sub>0</sub><sup>2</sup>β<sub>th</sub>ΔTHg / μ<sup>2</sup>L<sub>t</sub></td><td>0.90</td><td>0.90</td><td>2.2</td><td>(D<sup>29</sup>ρ<sub>10</sub><sup>2</sup>β<sub>th</sub>ΔTH<sup>8</sup>g / μ<sup>2</sup>L<sub>t</sub>)<sup>4/7</sup></td><td>0.80</td><td>0.80</td></tr></tbody></table>

## V.B. Residuals diagnosis

Although $Q_{CV}^2$ values are here close to 1, residual analysis provides an additional check of model validity under the assumed setting (first-degree polynomial with additive Gaussian noise). We focus here on the turbulent dataset with $\tau = 0.2$, corresponding to the most penalizing configuration. Figure 4(a) shows studentized residuals versus fitted values: they scatter randomly around $y = 0$, with no detectable pattern, indicating that the model $\hat{f}(x_{in}, \hat{\gamma}_{in})$ is well specified. This confirms not only the adequacy of the Gaussian noise assumption, but also the correctness of the chosen polynomial degree. Figure 4(b) presents a Q-Q plot against the Gaussian reference line; residuals align closely, confirming normality. These results are expected given the controlled test case, but they illustrate that the algorithm behaves consistently with its assumptions. Additional diagnostics, such as histograms of standardized residuals or formal normality tests (e.g. Shapiro-Wilk), could complement the graphical analysis. For practical applications, all these diagnostics remain essential. When the governing law is unknown and measurement errors affect the underlying physical variables, the Gaussian noise assumption and the choice of polynomial degree $P$ must be validated through residual analysis, to detect possible non-linearities or heteroscedasticity.

## V.C. Uncertainty intervals for the optimized parameters via bootstrap methods

Once the model assumptions have been validated, we now propose to evaluate confidence intervals (CIs) associated with the optimized parameters, in particular the exponents $\hat{\omega}_{in}$ (derived from $\hat{\gamma}_{in}$). Such CIs are essential to assess whether these parameters are truly identifiable from
<figure><img src="image_4.png" /></figure>

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

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

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

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

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

<figure><img src="image_10.png" /><figcaption>Fig. 3. Noisy simulated *Re* values (blue points) and predicted linear model (solid red line) plotted against the estimated dimensionless number, for different noise levels.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 4. Residual analysis for the turbulent dataset with $\tau = 0.2$.</figcaption></figure>

the available data, i.e. whether the dataset carries enough information to characterize them, and to account for residual uncertainty when interpreting the estimated dimensionless number. For the regression coefficient $\hat{\beta}(\hat{\gamma}_{in})$, CIs can be constructed directly under the Gaussian regression framework, yielding standard analytical bounds conditional on $\hat{\gamma}_{in}$. In contrast, no closed-form solution exists for CIs on parameters $\hat{\gamma}_{in}$, since they are obtained by cross-validation. Consequently, bootstrap methods provide a practical way to approximate the sampling distribution of $\hat{\gamma}_{in}$, and thereby of $\hat{\omega}_{in}$. CIs are thus obtained for $\hat{\omega}_{in}$ based on bootstrap replications of $\hat{\gamma}_{in}$, following the classical resampling approach introduced by Efron [32]:

1. Generate $B$ bootstrap samples $\{D_{N,b}^*\}_{b=1}^B$ from the original dataset $D_N$ using resampling with replacement. Each bootstrap sample has the same size as the original dataset.

2. Apply the D4NL method to each bootstrap sample to obtain $\hat{\gamma}_{in}^*(b)$. From these, derive the corresponding exponents $\hat{\omega}_{in}^*(b)$, yielding the collection $\{\hat{\omega}_{in}^*(b)\}_{b=1}^B$.

3. For each parameter $\omega_i$, compute summary statistics from the bootstrap replicates $\{\hat{\omega}_{in}^*(b)\}_{b=1}^B$, such as mean $\bar{\omega}_{in,i}^{(B)} = \frac{1}{B} \sum_{b=1}^B \hat{\omega}_{in,i}^*(b)$ and standard deviation $\hat{\sigma}_{in,i}^{(B)} = \left( \frac{1}{B-1} \sum_{b=1}^B \left[ \hat{\omega}_{in,i}^*(b) - \bar{\omega}_{in,i}^{(B)} \right]^2 \right)^{1/2}$.

4. Construct percentile bootstrap CI of level $(1 - \alpha)$ for each parameter $\omega_i$ by taking the empirical $\alpha/2$ and $1 - \alpha/2$ quantiles of the bootstrap distribution $\{\hat{\omega}_{in,i}^*(b)\}_{b=1}^B$. The resulting
bootstrap CI is:

$$ \left[ \hat{\omega}_{\text{in},i}^{*(\alpha/2)}, \hat{\omega}_{\text{in},i}^{*(1-\alpha/2)} \right] \qquad (32) $$

where $\hat{\omega}_{\text{in},i}^{*(p)}$ denotes the p-th empirical quantile of the bootstrap replicates.

In practice, the number of bootstrap replications B must be large enough to ensure stability of the estimates. A practical recommendation is to monitor the stabilization of the bootstrap variance $\hat{\sigma}_{\text{in},i}^{(B)}$ as B increases, and to adopt the smallest B for which changes become negligible. In this study, convergence was reached for B = 400, which was therefore used in all subsequent evaluations. The 95% bootstrap CIs obtained for the turbulent regime dataset with $\tau = 0.2$ are given in Table III. The estimated power and bootstrap samples are factorized by $\frac{4}{7}$ for identification purposes with the true physical model $Re = 2.87 \left( \frac{D^3 \rho_0^2 \beta_{\text{th}} \Delta T H g}{\mu^2 L_T} \right)^{4/7}$.

<table><caption>TABLE III
Bootstrap statistics for the turbulent regime dataset with $\tau = 0.2$.</caption><thead><tr><th>Index i</th><th>Associated variable</th><th>ω̂<sub>i</sub> in (ω̂<sub>i</sub>)<sup>4/7</sup></th><th>95% bootstrap CI</th><th>|σ̂<sup>(B)</sup><sub>in,i</sub>/ω̄<sup>(B)</sup><sub>in,i</sub>| in %</th></tr></thead><tbody><tr><td>1</td><td>D</td><td>2.97</td><td>[2.60; 3.38]</td><td>6.6</td></tr><tr><td>2</td><td>μ</td><td>-2.00</td><td>[-2.29; -1.74]</td><td>6.9</td></tr><tr><td>3</td><td>ρ<sub>0</sub></td><td>2.00</td><td>[1.74; 2.29]</td><td>6.9</td></tr><tr><td>4</td><td>β<sub>th</sub></td><td>1.00</td><td>[0.881; 1.12]</td><td>6.2</td></tr><tr><td>5</td><td>ΔT</td><td>1.00</td><td>[0.881; 1.12]</td><td>6.2</td></tr><tr><td>6</td><td>H</td><td>1.04</td><td>[0.787; 1.26]</td><td>11</td></tr><tr><td>7</td><td>L<sub>t</sub></td><td>-1.00</td><td>[-1.13; -0.874]</td><td>6.7</td></tr><tr><td>8</td><td>g</td><td>1.00</td><td>[0.869; 1.15]</td><td>6.9</td></tr></tbody></table>

When considering the absolute coefficient of variation of the bootstrapped values ($|\hat{\sigma}_{\text{in},i}^{(B)} / \bar{\omega}_{\text{in},i}^{(B)}|$), $\beta_{\text{th}}$ and $\Delta T$ are the most precisely identified, while the one linked to H exhibits the largest uncertainty. Uncertainty on estimated exponents reflect three main factors: (i) the influence of each input variable $X_i$ on the output $\Pi_{\text{out}}$, (ii) the informativeness of the dataset w.r.t this relationship, and (iii) the redundancy of physical dimensions across variables. The latter is directly linked to the dimension of $A_{\text{dim,in}}$ null space, derived from the dimensionless constraint. From Table III, it appears that exponents linked to dimensions shared by many inputs (e.g. length L) are harder to disentangle than those associated with unique dimensions (e.g. temperature $\Theta$).
## VI. EXTENSION OF THE D4NL METHOD TO UNKNOWN DIMENSIONLESS OUTPUT

The D4NL method can be extended to the case where the dimensionless output $\Pi_{out}$, previously known, is now unknown and must be evaluated. To our knowledge, such evolution has not been approached in other data-driven dimensionless learning studies. The objective of the extended D4NL is to identify both $\hat{\Pi}_{out}$ and $\hat{\Pi}_{in}$ and to estimate the polynomial governing law $\hat{f}$ from a dataset uniquely composed of $d+1$ physical variables $\mathbf{X}$. In this framework, both $\Pi_{out}$ and $\Pi_{in}$ are constructed from $\mathbf{X}$, except for a quantity of interest $X_{QoI}$ which is included only in $\Pi_{out}$. The assumptions for this extension of D4NL method are:

- The dimensionless output $\Pi_{out}$ is unknown.

- $\Pi_{out}$ depends on a known physical quantity of interest (QoI), $X_{QoI}$, which belongs to the input set $\mathbf{X}$ and must be preserved in the construction of $\Pi_{out}$.

- The dominant input dimensionless number $\hat{\Pi}_{in}$ is built from the remaining inputs and does not depend on $X_{QoI}$.

More formally, we define $\mathbf{X}_{out} \subseteq \mathbf{X}$ with $X_{QoI} \in \mathbf{X}_{out}$ and $\mathbf{X}_{in} \subseteq \mathbf{X} \setminus \{X_{QoI}\}$. Using the reparametrization built upon the kernel of the matrices $\mathbf{A}_{dim,out}$ and $\mathbf{A}_{dim,in}$ (i.e. the dimensional matrices associated with $\mathbf{X}_{out}$ and $\mathbf{X}_{in}$), the dimensionless numbers are then represented as $\Pi_{out}(\mathbf{X}_{out}, \gamma_{out})$ and $\Pi_{in}(\mathbf{X}_{in}, \gamma_{in})$ (cf. Eqs. (1) and (8) in Section II).

### VI.A. Nested optimization scheme to estimate ($\hat{\pi}_{out}$, $\hat{\pi}_{in}$)

To obtain expressions for $\hat{\pi}_{out}(\mathbf{x}_{out}, \hat{\gamma}_{out})$, $\hat{\pi}_{in}(\mathbf{x}_{in}, \hat{\gamma}_{in})$, and the associated $\hat{f}$, the loss to be minimized is:

$$ (\hat{\gamma}_{out}(\mathcal{D}_n), \hat{\gamma}_{in}(\mathcal{D}_n)) = \underset{\gamma_{out} \in \mathcal{D}_{\gamma_{out}}, \gamma_{in} \in \mathcal{D}_{\gamma_{in}}}{\operatorname{argmin}} \left( \mathcal{L}\left(\pi_{out}(\gamma_{out}), \hat{f}_{CV}(\gamma_{in})\right) \right), \quad (33) $$

where $\mathcal{D}_{\gamma_{out}}$ (resp. $\mathcal{D}_{\gamma_{in}}$) is the feasible domain of $\gamma_{out}$ (resp. $\gamma_{in}$) derived from hypercube constraints imposed on $\omega_{out}$ (resp. $\omega_{in}$) and projected on the kernel of $\mathbf{A}_{dim,out}$ (resp. $\mathbf{A}_{dim,in}$) (cf. Section III.C).
The solving strategy of Eq. (33) follows the same principles as in the original D4NL formulation with known $\Pi_{out}$, with an extension using a nested optimization scheme. In the outer loop, we optimize $\gamma_{out}$ within $\mathcal{D}_{\gamma_{out}}$; while in the inner loop, for each candidate $\gamma_{out}$, we optimize $\gamma_{in}$ within $\mathcal{D}_{\gamma_{in}}$ conditionally to $\gamma_{out}$. This corresponds to solving

$$ \hat{\gamma}_{\text{out}}(\mathcal{D}_n) = \underset{\gamma_{\text{out}} \in \mathcal{D}_{\gamma_{\text{out}}}}{\operatorname{argmin}} \left( \min_{\gamma_{\text{in}} \in \mathcal{D}_{\gamma_{\text{in}}}} \mathcal{L}(\pi_{\text{out}}(\gamma_{\text{out}}), \hat{f}_{\text{CV}}(\gamma_{\text{in}})) \right). \quad (34) $$

Once the optimal $\hat{\gamma}_{out}$ is identified, the corresponding inner optimum $\hat{\gamma}_{in}$ is obtained as:

$$ \hat{\gamma}_{\text{in}}(\mathcal{D}_n | \hat{\gamma}_{\text{out}}) = \underset{\gamma_{\text{in}} \in \mathcal{D}_{\gamma_{\text{in}}}}{\operatorname{argmin}} \mathcal{L}(\pi_{\text{out}}(\hat{\gamma}_{\text{out}}), \hat{f}_{\text{CV}}(\gamma_{\text{in}})). \quad (35) $$

This nested approach ensures that the optimal pair $(\hat{\gamma}_{out}, \hat{\gamma}_{in})$ is consistently identified, while preserving the robustness of the search against local minima thanks to the multistart initialization. Moreover, both loops rely on the same algorithmic structure as in standard D4NL: a multistart procedure is first applied to generate support points uniformly distributed over the feasible domain, and each starting point is then refined using a constrained descent algorithm to converge towards a local minimum. The nested scheme thus relies on two distinct sets of multistart support points: $N_{SP,out}$ for the outer loop optimization of $\gamma_{out}$, and $N_{SP,in}$ for the inner loop optimization of $\gamma_{in}$.

## VI.B. Additional ANOVA-based constraint for $\hat{\pi}_{out}$

To avoid pathological cases in which the outer loop inadvertently excludes $X_{QoI}$—yielding $\hat{\pi}_{out} = \hat{\pi}_{in}$ and a trivial $\hat{f} \equiv 1$—we propose to enforce an influence constraint on $X_{QoI}$ during the outer optimization. Concretely, we require that the contribution of $X_{QoI}$ to $\Pi_{out}$ is non-negligible, so that $X_{QoI}$ plays a significant role in the construction of $\hat{\Pi}_{out}$. This amounts to examine $\Pi_{out}$ functional decomposition as a function of its entries $X$, here considered as function parameters and not as random variables. We therefore only suppose the independence of the parameters (on their domain of definition $\mathcal{D}_X \in \mathbb{R}^{d+1}$) to analyze the contribution of $X_{QoI}$ to $\Pi_{out}$ function.

Under this assumption, the ANOVA (ANalysis Of VAriance) decomposition and associated Sobol' indices [33–35] provide a natural framework, since a candidate solution $\Pi_{out} := \prod_{i=1}^{d+1} X_i^{\omega_i}$ admits an explicit analytical form in terms of the $X_i$ (for fixed $\omega$), making analytical computation or Monte Carlo estimation of Sobol' indices straightforward. Importantly, we propose to consider
the ANOVA decomposition of $\log \Pi_{\text{out}}$ rather than of $\Pi_{\text{out}}$ itself. This choice is motivated by the multiplicative structure of $\Pi_{\text{out}} = \prod_{i=1}^{d+1} X_i^{\omega_i}$, which can distort the relative contributions of each $X_i$ when analyzed directly. Taking the logarithm linearizes the product into a sum so that the model becomes additive in the independent variables $\log X_i$. In this linear setting, the ANOVA decomposition is exact and the Sobol' indices admit closed-form expressions: the first-order and total indices coincide.

In practice, for each candidate vector $\gamma_{\text{out}}$ (and its associated exponents $\omega_{\text{out}}$), we can directly evaluate the total Sobol' index of $X_{\text{QoI}}$ on $\log \Pi_{\text{out}}$. Writing

$$ \log \Pi_{\text{out}} = \omega_{\text{out,QoI}} \log X_{\text{QoI}} + \sum_{\substack{i=1 \\ i \neq \text{QoI}}}^{d+1} \omega_{\text{out},i} \log X_i. \quad (36) $$

Since the function parameters $X_i$ ($i \in [1, d+1]$) are assumed to be independent of each other (for the analysis of functional dependence of $\Pi_{\text{out}}$ only), the variance decomposition yields:

$$ \mathrm{Var}(\log \Pi_{\text{out}}) = \omega_{\text{out,QoI}}^2 \mathrm{Var}(\log X_{\text{QoI}}) + \sum_{\substack{i=1 \\ i \neq \text{QoI}}}^{d+1} \omega_{\text{out},i}^2 \mathrm{Var}(\log X_i). \quad (37) $$

Denoting by $\mathbf{X}_{-\text{QoI}} := \mathbf{X} \setminus \{X_{\text{QoI}}\}$, the total Sobol' index of $X_{\text{QoI}}$, equal to the first-order index here, is then given analytically by:

$$ S_{\text{QoI}}^{\top}(\omega_{\text{out}}) = 1 - \frac{\mathrm{Var}(\mathbb{E}[\log \Pi_{\text{out}} | \mathbf{X}_{-\text{QoI}}])}{\mathrm{Var}(\log \Pi_{\text{out}})} \quad (38) $$

$$ = \frac{\mathrm{Var}(\mathbb{E}[\log \Pi_{\text{out}} | X_{\text{QoI}}])}{\mathrm{Var}(\log \Pi_{\text{out}})} \quad (39) $$

$$ = \frac{\omega_{\text{out,QoI}}^2 \mathrm{Var}(\log X_{\text{QoI}})}{\sum_i \omega_{\text{out},i}^2 \mathrm{Var}(\log X_i)}. \quad (40) $$

The feasible set for the outer optimization can then be restricted to:

$$ \mathcal{D}_{\gamma_{\text{out}}}^{\text{Sob}}(\eta) = \{\gamma_{\text{out}} \in \mathcal{D}_{\gamma_{\text{out}}} \mid S_{\text{QoI}}^{\top}(\omega_{\text{out}}) \geq \eta, \eta \in \mathbb{R}\}, \quad (41) $$

which enforces a minimal influence threshold $\eta$ for $X_{\text{QoI}}$ in $\hat{\Pi}_{\text{out}}$. The optimization problem in
Eq (34) becomes

$$ \hat{\gamma}_{\text{out}}(D_n, \eta) = \underset{\gamma_{\text{out}} \in \mathcal{D}_{\gamma_{\text{out}}}^{\text{Sob}}(\eta)}{\operatorname{argmin}} \left( \min_{\gamma_{\text{in}} \in \mathcal{D}_{\gamma_{\text{in}}}} \mathcal{L}\left(\pi_{\text{out}}(\gamma_{\text{out}}), \hat{f}_{\text{CV}}(\gamma_{\text{in}})\right) \right). \quad (42) $$

Note that this additional constraint is quadratic in $\gamma_{\text{out}}$, since the coefficients $\omega_{\text{out}}$ are obtained by projection $\omega_{\text{out}} = V_{\text{out}}\gamma_{\text{out}}$, and the Sobol' index reduces to a ratio of quadratic forms in these coefficients (see Appendix B). In practice, this inequality is added to the constraint set enforced by the constrained descent algorithm during loss minimization, for each multistart initialization.

## VI.C. Selection of solutions with high physical interpretability

To approach the canonical forms of dimensionless numbers commonly found in the literature, and thereby enhance the physical interpretability of the derived invariants, we introduce a post-processing step applied to candidate solutions. This step does not restrict the feasible domain during optimization, but rather acts afterwards on solutions that already satisfy the loss criterion and the feasibility constraints. The idea is to round the exponents of the dimensionless numbers to the nearest usual fractions, so that the resulting forms are closer to those traditionally used in physics and engineering practice.

Unlike the Sobol-based constraint, which restricts the optimization feasible prior domain, the rounding constraint is applied a posteriori to candidate solutions. Specifically, once solutions ($\hat{\gamma}_{\text{out}}$, $\hat{\gamma}_{\text{in}}$) have been identified from multistart procedure and satisfy the feasibility conditions, their associated exponent vectors $\hat{\omega}_{\text{out}} = V_{\text{out}}\hat{\gamma}_{\text{out}}$ and $\hat{\omega}_{\text{in}} = V_{\text{in}}\hat{\gamma}_{\text{in}}$ are rounded to the nearest usual fractions chosen among $\{\frac{1}{4}, \frac{1}{3}, \frac{1}{2}\}$ within $[-\omega_{\text{out}}^{\text{max}}, \omega_{\text{out}}^{\text{max}}]$ (with $\omega_{\text{out}}^{\text{max}} = \max(\omega_{\text{out}}^{\text{max}}, \omega_{\text{in}}^{\text{max}}) \in \mathbb{R}$). The resulting rounded vectors are denoted $\hat{\omega}_{\text{out}}^{\text{uf}}$ and $\hat{\omega}_{\text{in}}^{\text{uf}}$ (with “uf” standing for usual fractions). Only solutions that remain dimensionless after rounding are retained in the final admissible set:

$$ S^{\text{uf}} = \{(\hat{\gamma}_{\text{out}}, \hat{\gamma}_{\text{in}}) \mid \mathbf{A}_{\text{dim,out}}\hat{\omega}_{\text{out}}^{\text{uf}} = \mathbf{0}, \mathbf{A}_{\text{dim,in}}\hat{\omega}_{\text{in}}^{\text{uf}} = \mathbf{0}\}. \quad (43) $$

The D4NL estimation and validation procedures extended to $\Pi_{\text{out}}$ unknown are summarized by Algorithms 3 and 4 in Appendix A.
VI.D. Application to the thermosyphon test case

The extended D4NL method is now applied to the single-phase natural circulation (thermosyphon) test case in the laminar regime. The goal is to recover the laminar form of Eq. (29), expressed as:

$$ \frac{\rho_0 v D}{\mu} = \beta_{\text{true}} \frac{D^3 \rho_0^2 \beta_{\text{th}} \Delta T H g}{\mu^2 L_t}, \qquad (44) $$

or equivalently, to identify governing relation expressed in terms of dimensionless numbers analogous to the Reynolds number ($Re = \frac{\rho_0 v D}{\mu}$) and the modified Grashof number ($Gr_m = \frac{D^3 \rho_0^2 \beta_{\text{th}} \Delta T H g}{\mu^2 L_t}$).

In addition to the set of physical variables previously considered, the velocity $v$ (computed from Eq. (22)) is explicitly included, so that $X = (D, \mu, \rho_0, \beta_{\text{th}}, \Delta T, H, L_t, g, v)^{\top}$. The outer dimensionless group $\Pi_{\text{out}}$ is constructed from the full set $X_{\text{out}} = X$, with $v$ designated as the quantity of interest (QoI) to be preserved. In contrast, the inner group $\Pi_{\text{in}}$ is defined on the reduced set $X_{\text{in}} = X \setminus \{v\}$, excluding the velocity. Because $v$ is itself a deterministic function of $X_{\text{in}}$ (via Eq. (22)), this configuration may introduce identifiability issues in the optimization. However, the Sobol-based constraint on the functional dependence of $\Pi_{\text{out}}$ on $X_{\text{QoI}}$ acts as a regularization mechanism, thereby mitigating the compensation effects among dependent variables and stabilizing the recovered exponents.

As in Section IV.B, we generate a dataset of $N = 200$ i.i.d. samples $\{x^{(j)}\}_{j=1}^{200}$. Unlike the previous case, the datasets used here are noiseless, which allows us to focus on the structural recovery of governing laws without the confounding effect of measurement uncertainty. Regarding the dimensional structure, the fundamental dimension matrix associated with $\Pi_{\text{in}}$ remains unchanged and is equal to $A_{\text{dim,in}}$ given in Equation (31). By contrast, the matrix associated with $\Pi_{\text{out}}$, which now includes the velocity $v$, becomes

$$ A_{\text{dim,out}} = \begin{pmatrix} D & \mu & \rho_0 & \beta_{\text{th}} & \Delta T & H & L_t & g & v \\ 1 & -1 & -3 & 0 & 0 & 1 & 1 & 1 & 1 \\ 0 & -1 & 0 & 0 & 0 & 0 & 0 & -2 & -1 \\ 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 & 1 & 0 & 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} L \\ T \\ M \\ \Theta \end{pmatrix}. \qquad (45) $$
$A_{\text{dim,in}} \in \mathbb{R}^{4 \times 8}$ when $A_{\text{dim,out}} \in \mathbb{R}^{4 \times 9}$.

To perform the nested optimization, we allocate more support points to the multistart procedure in the outer loop than in the inner one. From the application with known $\Pi_{\text{out}}$, we observed that only a limited number of support points is required for the inner optimization. We therefore concentrate the exploration effort on the outer parameter space $D_{\gamma_{\text{out}}}$, where the algorithmic behavior is less predictable. In practice, we set $N_{\text{SP,out}} \approx 10 N_{\text{SP,in}}$. The same convergence tolerance is used in both loops. Table IV reports the solutions in $S^{\text{uf}}$ obtained with $N_{\text{SP,out}} = 200$.

<table><caption>TABLE IV
Solutions obtained with the nested optimization scheme.</caption><thead><tr><th>Index</th><th>Equation</th><th>Simplified equation</th><th>Loss value</th></tr></thead><tbody><tr><td>1</td><td>$\frac{L_t^{\frac{3}{4}} v^{\frac{3}{2}}}{\sqrt{D} H \sqrt{\beta_{\text{th}} g^{\frac{3}{4}} \sqrt{\Delta T}}} = \hat{\beta}_1 \frac{D^{\frac{5}{2}} \sqrt{H} \beta_{\text{th}} g^{\frac{3}{4}} \rho_0^{\frac{3}{2}} \Delta T}{L_t^{\frac{3}{4}} \mu^{\frac{3}{2}}}$</td><td>$\frac{L_t^{\frac{3}{2}} \mu^{\frac{3}{2}} v^{\frac{3}{2}}}{D^3 H^{\frac{3}{2}} \beta_{\text{th}} g^{\frac{3}{2}} \rho_0^{\frac{3}{2}} \Delta T^{\frac{3}{2}}} = \hat{\beta}_1$</td><td>$3.5 \times 10^{-19}$</td></tr><tr><td>2</td><td>$\frac{H^{\frac{2}{3}} \sqrt{L_t} \beta_{\text{th}}^{\frac{5}{3}} \sqrt{\rho_0 v^{\frac{3}{2}}} \Delta T^{\frac{5}{4}}}{D^{\frac{3}{4}} \sqrt{g} \sqrt{\mu}} = \hat{\beta}_2 \frac{D^{\frac{9}{4}} H^{\frac{9}{4}} \beta_{\text{th}}^{\frac{11}{4}} g \rho_0^{\frac{2}{3}} \Delta T^{\frac{11}{4}}}{L_t^{\frac{3}{4}} \mu^2}$</td><td>$\frac{L_t^{\frac{19}{12}} \mu^{\frac{3}{2}} v^{\frac{3}{2}}}{D^3 H^{\frac{19}{12}} \beta_{\text{th}} g^{\frac{3}{2}} \rho_0^{\frac{3}{2}} \Delta T^{\frac{3}{2}}} = \hat{\beta}_2$</td><td>$4.9 \times 10^{-4}$</td></tr><tr><td>3</td><td>$\frac{\sqrt[4]{H} \sqrt[4]{v}}{\sqrt[4]{D} \sqrt[4]{L_t} \sqrt[4]{g}} = \hat{\beta}_3 \frac{\sqrt[4]{D} \sqrt[4]{H} \sqrt[4]{\beta_{\text{th}}} \sqrt[4]{\rho_0} \sqrt[4]{\Delta T}}{\sqrt[4]{L_t} \sqrt[4]{\mu}}$</td><td>$\frac{\sqrt[4]{L_t} \sqrt[4]{\mu} \sqrt[4]{v}}{\sqrt{D} \sqrt[4]{H} \sqrt[4]{\beta_{\text{th}}} \sqrt[4]{g} \sqrt[4]{\rho_0} \sqrt[4]{\Delta T}} = \hat{\beta}_3$</td><td>$1.5 \times 10^{-3}$</td></tr><tr><td>4</td><td>$\frac{H^{\frac{4}{3}} v}{D^{\frac{3}{4}} L_t^{\frac{3}{4}} \sqrt[4]{\beta_{\text{th}}} \sqrt{g} \sqrt[4]{\Delta T}} = \hat{\beta}_4 \frac{D^{\frac{5}{4}} H^{\frac{7}{3}} \beta_{\text{th}}^{\frac{3}{4}} \sqrt{g} \rho_0 \Delta T^{\frac{3}{4}}}{L_t^{\frac{7}{4}} \mu}$</td><td>$\frac{L_t \mu v}{D^2 H \beta_{\text{th}} g \rho_0 \Delta T} = \hat{\beta}_4$</td><td>$3.3 \times 10^{-19}$</td></tr></tbody></table>

To assess whether the governing laws in the second column of Table IV are similar with Equation (44), we compute their simplified forms (with redundant variables removed) in the third column using the SymPy package [36]. The simplification process applied to Eq. (44) yields

$$ \frac{L_t \mu v}{D^2 H \beta_{\text{th}} g \rho_0 \Delta T} = \beta_{\text{true}}. \qquad (46) $$

which coincides exactly with simplified equation number 4 in Table IV, matches equation 1 up to a power 2/3. Both solutions 1 and 4 yield very low loss values: the algorithm successfully recovered governing laws which are similar to Eq. (44) with new $\hat{\pi}_{\text{in}}$ and $\hat{\pi}_{\text{out}}$.

## VII. CONCLUSION

The licensing of nuclear passive safety systems requires scaling analyses to support the design of reduced size experiments that reproduce the relevant safety related phenomena. However, scale distortions are unavoidable in such experiments, and their impact on the predicted system
performance is often difficult to quantify. These distortions can nevertheless be mitigated during
the design phase if the loop is built around a dominant dimensionless number. This strategy is
commonly applied when the governing dimensionless numbers are known, such as the Reynolds
number in thermal hydraulics. Recently, data driven methods for learning dimensionless num-
bers have emerged, offering new perspectives for designing test facilities by identifying dominant
dimensionless groups directly from available data.

To begin with, this paper presents a new dimensionless-number learning method, named D4NL for Data-Driven Discovery of Dimensionless Numbers and governing Laws. This polynomial regression approach benefits from a clear mathematical formulation combined with sophisticated machine learning and optimization tools. Applied to a simulated single phase natural circulation dataset, the method successfully recovers the governing laws that link the known Reynolds number to the yet-to-be-identified modified Grashof number, in both laminar and turbulent regimes. As a perspective, laminar and turbulent flow regimes could be treated in a single dataset, leveraging future algorithm capabilities in regime detection.

Mostly, the D4NL formalization paves the way for extending the method to multivariate
and nonlinear governing-law identification. The generalized approach will be tested on aggregated
experimental data involving nonlinear phenomena, like two-phase flow. Evaluating the method on
experimental data will enable to validate the structural modeling choices, including the assumed
error distribution, and to assess the algorithm's efficiency w.r.t. the ratio of observations to
variables.

In addition, this paper presents the extension of the D4NL method to the case where the
output dimensionless number is unknown. In this setting, the user only needs to specify the
quantity of interest that must be preserved in the output dimensionless number. The algorithm
relies on a nested optimization procedure with additional constraints designed to ease convergence
towards relevant dimensionless output and input. This extended version is tested on single-phase
natural circulation data and identifies governing laws based on variants of Reynolds and Grashof
numbers. The results highlight the method's strong potential for future use across a wide range of
physical phenomena, relevant to the design of engineering system test facilities.
A. ALGORITHMS OF THE D4NL METHOD

This appendix provides the detailed algorithms used in the D4NL method. For clarity, it is organized into four parts: (1) estimation of $\gamma_{in}$ via K-Fold CV when $\Pi_{out}$ is known, (2) validation by nested cross-validation in the known $\Pi_{out}$ case, (3) estimation of $\gamma_{out}$ and $\gamma_{in}$ when $\Pi_{out}$ is unknown, with Sobol and rounding constraints, and (4) validation by outer cross-validation in the unknown $\Pi_{out}$ case.

**Algorithm 1** Estimation procedure of D4NL when $\Pi_{out}$ is known: estimation of $\gamma_{in}$ via K-Fold Cross-Validation

**Require:** Dataset $\mathcal{D}_n = \{(\mathbf{x}_{in}^{(j)}, \pi_{out}^{(j)})\}_{j=1}^n$, K folds for CV, $N_{SP,in}$ support points for $\gamma_{in}$.

Compute the dimensionless basis $\mathbf{V}_{in}$ associated with $\mathbf{X}_{in}$.

Generate $N_{SP,in}$ starting points for $\gamma_{in}$ within the domain $\mathcal{D}_{\gamma_{in}}$.

Construct $\mathcal{L}_{CV}$, the loss function based on K-Fold CV (defined in Eqs. (13)–(11)).

for $i \leftarrow 1$ to $N_{SP,in}$ do

Initialize $\gamma_{in,i}$ at the $i$-th support point.

Run constrained gradient-based optimization from $\gamma_{in,i}$ to minimize $\mathcal{L}_{CV}$.

Save the final $\hat{\gamma}_{in,i}$ and its loss value.

end for

Select $\hat{\gamma}_{in}^*$ with the lowest $\mathcal{L}_{CV}$.

Compute $\hat{\omega}_{in}^* = \mathbf{V}_{in}\hat{\gamma}_{in}^*$.

Form $\hat{\pi}_{in}(\mathbf{x}_{in}, \hat{\omega}_{in}^*)$ and compute $\hat{\beta}(\hat{\gamma}_{in}^*)$ on $\mathcal{D}_n$.

return $(\hat{\gamma}_{in}^*, \hat{\omega}_{in}^*, \hat{\pi}_{in}, \hat{\beta}(\hat{\gamma}_{in}^*))$.

**Algorithm 2** Validation procedure of D4NL when $\Pi_{out}$ is known: nested cross-validation and final model construction

**Require:** Dataset $\mathcal{D}_n = \{(\mathbf{x}_{in}^{(j)}, \pi_{out}^{(j)})\}_{j=1}^n$, K' folds for outer CV, K folds for inner CV, $N_{SP,in}$ support points for $\gamma_{in}$.

for $k' \leftarrow 1$ to $K'$ do

Define outer training set $\mathcal{D}_{n\setminus\{k'\}}$.

Call Algorithm 1 on dataset $\mathcal{D}_{n\setminus\{k'\}}$ with K folds and $N_{SP,in}$ support points.

Obtain $(\hat{\gamma}_{in,k'}, \hat{\omega}_{in,k'}, \hat{\pi}_{in,k'}, \hat{\beta}_{k'})$.

Predict outputs on the left-out fold $k'$ and save residuals $\pi_{out}^{(i)} - \hat{\pi}_{out}^{(i)}$ for all $j \in k'$.

end for

Aggregate residuals across all $K'$ folds and deduce $Q_{CV}^2$ via Eq. (17).

**Final model construction:**

Call Algorithm 1 on the full dataset $\mathcal{D}_n$ with K folds and $N_{SP,in}$ support points.

Obtain $(\hat{\gamma}_{in}^*, \hat{\omega}_{in}^*, \hat{\pi}_{in}, \hat{\beta}^*)$.

Define the final governing law estimator $\hat{f}(\mathbf{x}_{in}, \hat{\gamma}_{in}^*)$.

return $(\hat{\gamma}_{in}^*, \hat{\omega}_{in}^*, \hat{\pi}_{in}, \hat{f}(\cdot), Q_{CV}^2)$.
**Algorithm 3** Estimation procedure of D4NL when $\Pi_{\text{out}}$ is unknown: nested optimization of $\gamma_{\text{out}}$ and $\gamma_{\text{in}}$ with Sobol and rounding constraints

**Require:** Dataset $\mathcal{D}_n = \{\mathbf{x}^{(j)}\}_{j=1}^n$, quantity of interest $X_{\text{QoI}}$, K folds for inner CV, $N_{\text{SP,in}}$ support points for $\gamma_{\text{in}}$, $N_{\text{SP,out}}$ support points for $\gamma_{\text{out}}$, Sobol threshold $\eta$.

Define $\mathbf{X}_{\text{out}} \subseteq \mathbf{X}$ with $X_{\text{QoI}} \in \mathbf{X}_{\text{out}}$.

Define $\mathbf{X}_{\text{in}} \subseteq \mathbf{X} \setminus \{X_{\text{QoI}}\}$.

Compute the dimensionless basis $\mathbf{V}_{\text{out}}$ associated with $\mathbf{A}_{\text{dim,out}}$ of $\mathbf{X}_{\text{out}}$.

Compute the dimensionless basis $\mathbf{V}_{\text{in}}$ associated with $\mathbf{A}_{\text{dim,in}}$ of $\mathbf{X}_{\text{in}}$.

Generate $N_{\text{SP,out}}$ starting points for $\gamma_{\text{out}}$ within the domain $\mathcal{D}_{\gamma_{\text{out}}}$.

**for** i ← 1 to $N_{\text{SP,out}}$ **do**

Initialize candidate $\gamma_{\text{out},i}$.

Run the outer constrained optimization as in Eq. (42), enforcing $S_{\text{QoI}}^{\top}(\omega_{\text{out}}) \ge \eta$ with $\omega_{\text{out}} = \mathbf{V}_{\text{out}}\gamma_{\text{out}}$: within each outer iteration, call Algorithm 1 to solve the inner minimization on $\gamma_{\text{in}}$ conditionally to the current $\gamma_{\text{out}}$, using K, $N_{\text{SP,in}}$, the dataset restricted to $\mathbf{X}_{\text{in}}$, and the corresponding $\pi_{\text{out}}(\gamma_{\text{out}})$.

Recover $\hat{\gamma}_{\text{in}}(\hat{\gamma}_{\text{out},i})$ and compute the associated $\hat{\beta}^{(i)}$.

Save $(\hat{\gamma}_{\text{out},i}, \hat{\gamma}_{\text{in}}(\hat{\gamma}_{\text{out},i}), \hat{\beta}^{(i)})$ and its loss value.

**end for**

Select $\hat{\gamma}_{\text{out}}^*$ that minimizes the loss among admissible rounded candidates in $\mathcal{S}^{\text{uf}}$ (cf. Eq (43)).

Recover $\hat{\gamma}_{\text{in}}^* = \hat{\gamma}_{\text{in}}(\hat{\gamma}_{\text{out}}^*)$ and associated $\hat{\beta}^*$.

Define the final governing law estimator $\hat{f}(\mathbf{x}; \hat{\gamma}_{\text{out}}^*, \hat{\gamma}_{\text{in}}^*, \hat{\omega}_{\text{out}}^{\text{uf}}, \hat{\omega}_{\text{in}}^{\text{uf}})$.

**return** $(\hat{\gamma}_{\text{out}}^*, \hat{\gamma}_{\text{in}}^*, \hat{\omega}_{\text{out}}^{\text{uf}}, \hat{\omega}_{\text{in}}^{\text{uf}}, \hat{\pi}_{\text{out}}, \hat{\pi}_{\text{in}}, \hat{\beta}^*, \hat{f}(\cdot))$.

**Algorithm 4** Validation procedure of the D4NL when $\Pi_{\text{out}}$ is unknown: outer cross-validation and final model construction

**Require:** Dataset $\mathcal{D}_n = \{\mathbf{x}^{(j)}\}_{j=1}^n$, K' folds for outer CV, Sobol threshold $\eta$, parameters $N_{\text{SP,in}}$, $N_{\text{SP,out}}$, K.

**for** each fold k = 1 to K' **do**

Split $\mathcal{D}_n$ into training set $\mathcal{D}_n^{(-k)}$ and validation set $\mathcal{D}_n^{(k)}$.

Apply Algorithm 3 on $\mathcal{D}_n^{(-k)}$ with inputs $(K, N_{\text{SP,in}}, N_{\text{SP,out}}, \eta)$.

Obtain estimated parameters $(\hat{\gamma}_{\text{out}}^{(-k)}, \hat{\gamma}_{\text{in}}^{(-k)}, \hat{\beta}^{(-k)}, \hat{f}^{(-k)})$.

Evaluate $\hat{f}^{(-k)}$ on the validation set $\mathcal{D}_n^{(k)}$ and compute the validation loss.

Save the validation loss for fold k.

**end for**

Aggregate the K' validation losses to obtain the final performance metrics (mean, variance, distribution).

**return** aggregated validation results and the final estimator trained on the full dataset $\mathcal{D}_n$.

**return** $(\hat{\gamma}_{\text{out}}^*, \hat{\gamma}_{\text{in}}^*, \hat{\omega}_{\text{out}}^{\text{uf}}, \hat{\omega}_{\text{in}}^{\text{uf}}, \hat{\pi}_{\text{out}}, \hat{\pi}_{\text{in}}, \hat{\beta}^*, \hat{f}(\cdot)$, validation metrics).
## B. REFORMULATION OF SOBOL-BASED INEQUALITY CONSTRAINT

Recalling that the Sobol index of the QoI is

$$S_{\text{QoI}}^{\top}(\omega_{\text{out}}) = \frac{\omega_{\text{out, QoI}}^2 \text{Var}(\log X_{\text{QoI}})}{\sum_{i} \omega_{\text{out, i}}^2 \text{Var}(\log X_i)},$$

and that $\omega_{\text{out}}$ is obtained from $\gamma_{\text{out}}$ via the projection

$$\omega_{\text{out}} = V_{\text{out}}\gamma_{\text{out}},$$

the inequality constraint

$$S_{\text{QoI}}^{\top}(\omega_{\text{out}}) \geq \eta$$

can be rewritten explicitly as

$$(\mathbf{v}_{\text{QoI}}^{\top}\gamma_{\text{out}})^2 \text{Var}(\log X_{\text{QoI}}) \geq \eta \sum_{i} (\mathbf{v}_{i}^{\top}\gamma_{\text{out}})^2 \text{Var}(\log X_i), \quad (47)$$

where $\mathbf{v}_i$ denotes the row of $V_{\text{out}}$ associated with variable $X_i$.

Equation (47) is thus a quadratic inequality in $\gamma_{\text{out}}$.
REFERENCES

[1] A. N. NAHAVANDI, F. S. CASTELLANA, and E. N. MORADKHANIAN, "Scaling Laws for Modeling Nuclear Reactor Systems," *Nuclear Science and Engineering*, **72**, 1, 75-83 (1979); 10.13182/nse79-a19310., URL http://dx.doi.org/10.13182/nse79-a19310.

[2] A. CIECHOCKI, S. CARNEVALI, D. BESTION, and L. ROSSI, "Scaling analyses of a SB-LOCA counterpart test between BETHSY and LSTF facilities and a three loops PWR," *Nuclear Engineering and Design*, **400**, 112051 (2022); 10.1016/j.nucengdes.2022.112051., URL http://dx.doi.org/10.1016/j.nucengdes.2022.112051.

[3] M. ISHII and I. KATAOKA, *Similarity analysis and scaling criteria for LWRs under single-phase and two-phase natural circulation* (1983); 10.2172/6312011., URL http://dx.doi.org/10.2172/6312011.

[4] M. ISHII, S. RAVANKAR, and R. DOWLATI, *Scientific design of Purdue University Multi-Dimensional Integral Test Assembly (PUMA) for GE SBWR* (1996); 10.2172/220555., URL http://dx.doi.org/10.2172/220555.

[5] M. ISHII, S. REVANKAR, T. LEONARDI, R. DOWLATI, M. BERTODANO, I. BABELLI, W. WANG, H. POKHARNA, V. RANSOM, R. VISKANTA, and J. HAN, "The three-level scaling approach with application to the Purdue University Multi-Dimensional Integral Test Assembly (PUMA)," *Nuclear Engineering and Design*, **186**, 1-2, 177-211 (1998); 10.1016/S0029-5493(98)00222-2., URL http://dx.doi.org/10.1016/S0029-5493(98)00222-2.

[6] N. ZUBER, "A Hierarchical, Two-tiered Scaling Analysis, Appendix D of An Integrated Structure and Scaling Methodology for Severe Accident Technical Issue Resolution," NUREG/CR 5809, U.S. Nuclear Regulatory Commission (1991).

[7] N. ZUBER, G. E. WILSON, M. ISHII, W. WULFF, B. BOYACK, A. DUKLER, P. GRIFFITH, J. HEALZER, R. HENRY, J. LEHNER, S. LEVY, F. MOODY, M. PILCH, B. SEHGAL, B. SPENCER, T. THEOFANOUS, and J. VALENTE, "An integrated structure and scaling methodology for severe accident technical issue resolution: Development of methodology," *Nuclear Engineering and Design*, **186**, 1-2, 1-21 (1998); 10.1016/S0029-5493(98)00215-5., URL http://dx.doi.org/10.1016/S0029-5493(98)00215-5.
[8] P. S. LYKOUDIS, "Non-dimensional numbers as ratios of characteristic times," *International Journal of Heat and Mass Transfer*, **33**, 7, 1568 (1990); 10.1016/0017-9310(90)90053-W.

[9] E. BUCKINGHAM, "The Principle of Similitude," *Nature*, **96**, 2406, 396 (1915); 10.1038/096396d0.

[10] P. F. MENDEZ and F. ORDÓÑEZ, "Scaling Laws From Statistical Data and Dimensional Analysis," *Journal of Applied Mechanics*, **72**, 5, 648 (2005); 10.1115/1.1943434.

[11] P. G. CONSTANTINE, Z. DEL ROSARIO, and G. IACCARINO, "Data-driven dimensional analysis: algorithms for unique and relevant dimensionless groups," (2017)URL https://arxiv.org/pdf/1708.04303.

[12] L. JOFRE, Z. R. DEL ROSARIO, and G. IACCARINO, "Data-driven dimensional analysis of heat transfer in irradiated particle-laden turbulent flow," *International Journal of Multiphase Flow*, **125**, 103198 (2020); 10.1016/j.ijmultiphaseflow.2019.103198.

[13] J. BAKARJI, J. CALLAHAM, S. L. BRUNTON, and J. N. KUTZ, "Dimensionally Consistent Learning with Buckingham Pi," (2022)URL https://arxiv.org/pdf/2202.04643.

[14] W. TENACHI, R. IBATA, and F. I. DIAKOGIANNIS, "Deep Symbolic Regression for Physics Guided by Units Constraints: Toward the Automated Discovery of Physical Laws," *The Astrophysical Journal*, **959**, 2, 99 (2023); 10.3847/1538-4357/ad014c.

[15] X. XIE, A. SAMAEI, J. GUO, W. K. LIU, and Z. GAN, "Data-driven discovery of dimensionless numbers and governing laws from scarce measurements," URL https://www.nature.com/articles/s41467-022-35084-w.

[16] W. FORD, "Gram-Schmidt Orthonormalization," *Numerical Linear Algebra with Applications*, 281-297, Elsevier; 10.1016/B978-0-12-394435-1.00014-4.

[17] T. HASTIE, R. TIBSHIRANI, and J. H. FRIEDMAN, *The elements of statistical learning: data mining, inference, and prediction*, Springer, New York (2009).

[18] E. FEKHARI, B. IOOSS, J. MURÉ, L. PRONZATO, and J. RENDAS, "Model predictivity assessment: incremental test-set selection and accuracy evaluation," N. SALVATI, C. PERNA,
S. MARCHETTI, and R. CHAMBERS (Editors), *Studies in Theoretical and Applied Statistics, SIS 2021, Pisa, Italy, June 21-25*, 315–347, Springer (2023).

[19] S. MAK and V. R. JOSEPH, “Support points,” *The Annals of Statistics*, **46**, 6A, 2562 (2018); 10.1214/17-AOS1629.

[20] C. DEMAY, B. IOOSS, L. LE GRATIET, and A. MARREL, “Model selection based on validation criteria for Gaussian process regression: An application with highlights on the predictive variance,” *Quality and Reliability Engineering International*, **38**, 3, 1482 (2022); 10.1002/qre.2973.

[21] V. KOUHIA, V. RIIKONEN, O.-P. KAUPPINEN, J. TELKKÄ, and J. HYVÄRINEN, “PASI – A test facility for research on passive heat removal,” *Nuclear Engineering and Design*, **383**, 111417 (2021); 10.1016/j.nucengdes.2021.111417.

[22] A. MANERA and T. H. J. J. VAN DER HAGEN, “Stability of Natural-Circulation-Cooled Boiling Water Reactors during Startup: Experimental Results,” *Nuclear Technology*, **143**, 1, 77 (2003); 10.13182/NT03-A3399.

[23] S. RENAUDIÈRE DE VAUX, P. AUBERT, B. GROSJEAN, and L. ROSSI, “Flow reversals in a natural circulation loop at atmospheric pressure,” *International Journal of Heat and Mass Transfer*, **235**, 126119 (2024); 10.1016/j.ijheatmasstransfer.2024.126119.

[24] P. VIJAYAN, M. SHARMA, and D. SAHA, “Steady state and stability characteristics of single-phase natural circulation in a rectangular loop with different heater and cooler orientations,” *Experimental Thermal and Fluid Science*, **31**, 8, 925 (2007); 10.1016/j.expthermflusci.2006.10.003.

[25] C. TOMPKINS and M. CORRADINI, “Flow pattern transition instabilities in a natural circulation cooling facility,” *Nuclear Engineering and Design*, **332**, 267 (2018); 10.1016/j.nucengdes.2018.03.018.

[26] S. JUNG, I. JENTZ, W. JUNG, B. LINDLEY, M. ANDERSON, and M. CORRADINI, “Reactor Vessel Cooling System for the GA-EMS Fast Modular Reactor: Scaled Experiments and Analysis,” *Advances in Thermal Hydraulics*, ATH 2024, Orlando, FL, 614–623 (2024).
[27] P. DENÉ, M. MONTOUT, E. GARCIA, J. TELKKÄ, V. RIIKONEN, and F. DAVID, "Frequency and amplitude of flashing-induced instability in an open natural circulation loop," *Nuclear Engineering and Design*, **424**, 113279 (2024); 10.1016/j.nucengdes.2024.113279.

[28] T. ZHANG and C. S. BROOKS, "Stability Tests and Analysis of a Low-Pressure Natural Circulation Loop with Flashing Instability," *Nuclear Technology*, **209**, 10, 1414 (2023); 10.1080/00295450.2022.2151823.

[29] P. K. VIJAYAN, A. K. NAYAK, and N. KUMAR, *Single-phase, two-phase and supercritical natural circulation systems*, Woodhead Publishing series in energy, Woodhead Publishing, an imprint of Elsevier, Duxford Cambridge, MA Kidlington (2019).

[30] J. BOUSSINESQ, *Essai sur la théorie des eaux courantes*, Gauthier-Villars (1897).

[31] X. FANG, Y. XU, and Z. ZHOU, "New correlations of single-phase friction factor for turbulent pipe flow and evaluation of existing single-phase friction factor correlations," *Nuclear Engineering and Design*, **241**, 3, 897 (2011); 10.1016/j.nucengdes.2010.12.019.

[32] B. EFRON and R. TIBSHIRANI, *An Introduction to the Bootstrap*, Chapman and Hall (1994); 10.1201/9780429246593.

[33] I. SOBOL', "Global sensitivity indices for non linear mathematical models and their Monte Carlo estimates," *Mathematics and Computers in Simulation*, **55**, 271 (2001).

[34] C. PRIEUR and S. TARANTOLA, "Variance-Based Sensitivity Analysis: Theory and Estimation Algorithms," R. GHANEM, D. HIGDON, and H. OWHADI (Editors), *Handbook of Uncertainty Quantification*, 1-23, Springer International Publishing.

[35] S. DA VEIGA, F. GAMBOA, B. IOOSS, and C. PRIEUR, *Basics and Trends in Sensitivity Analysis*, Society for Industrial and Applied Mathematics, Philadelphia, PA (2021); 10.1137/1.9781611976694.

[36] A. MEURER, C. P. SMITH, M. PAPROCKI, O. ČERTÍK, S. B. KIRPICHEV, M. ROCKLIN, A. KUMAR, S. IVANOV, J. K. MOORE, S. SINGH, T. RATHNAYAKE, S. VIG, B. E. GRANGER, R. P. MULLER, F. BONAZZI, H. GUPTA, S. VATS, F. JOHANSSON, F. PEDREGOSA, M. J. CURRY, A. R. TERREL, V. ROUČKA, A. SABOO, I. FERNANDO, S. KULAL, R. CIMRMAN,
and A. SCOPATZ, “SymPy: symbolic computing in Python,” *PeerJ Computer Science*, **3**, 103 (2017); 10.7717/peerj-cs.103.