G Geraci
email: ggeraci@stanford.edu

P M Congedo
email: pietro.congedo@inria.fr

R Abgrall

G Iaccarino

High-order statistics in global sensitivity analysis: decomposition and model reduction

Keywords: sensitivity analysis, high-order statistics, skewness, kurtosis, polynomial chaos, model reduction, ORCs turbines

ANalysis Of VAriance (ANOVA) is a common technique for computing a ranking of the input parameters in terms of their contribution to the output variance. Nevertheless, the variance is not an universal criterion for ranking variables, since non symmetric outputs could require higher order statistics for their description and analysis. In this work, we illustrate how third and fourth-order moments, i.e. skewness and kurtosis, respectively, can be decomposed mimicking the ANOVA approach. It is also shown how this decomposition is correlated to a Polynomial Chaos (PC) expansion leading to a simple strategy to compute each term. New sensitivity indices, based on the contribution to the skewness and kurtosis, are proposed. The outcome of the proposed analysis is depicted by considering several test functions. Moreover, the ranking of the sensitivity indices is shown to vary according to their statistics order. Furthermore, the problem of formulating a truncated polynomial representation of the original function is treated. Both the reduction of the number of dimensions and the reduction of the order of interaction between parameters are considered. In both cases, the impact on the reduction is assessed in terms of statistics, namely the probability density function. Feasibility of the proposed analysis in a real-case is then demonstrated by presenting the sensitivity analysis of the performances of a turbine cascade in an Organic Rankine Cycles (ORCs), in the presence of complex thermodynamic models and multiple sources of uncertainty.

Introduction

Optimization and design in the presence of uncertain operating conditions, material properties and manufacturing tolerances poses a tremendous challenge to the scientific computing community. Uncertainty quantification (UQ) approaches represent the inputs as random variables and seek to construct a statistical characterization of the quantities of interest. Several methodologies are proposed to tackle this issue, among those stochastic spectral methods [START_REF] Xiu | The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations[END_REF][START_REF] Xiu | High-Order Collocation Methods for Differential Equations with Random Inputs[END_REF][START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF][START_REF] Ghanem | Stochastic Finite Elements. A spectral approach[END_REF][START_REF] Nobile | A sparse grid stochastic collocation method for partial differential equations with random input data[END_REF], that can provide considerable speed-up in computational time when compared to Monte Carlo (MC) simulation. The presence of a large number of uncertain inputs leads to an exponential increase of the cost thus making these methodologies unfeasible [START_REF] Foo | Multi-element probabilistic collocation method in high dimensions[END_REF]. This situation becomes even more challenging when robust design optimization is tackled [START_REF] Kim | Efficient Shape Optimization Under Uncertainty Using Polynomial Chaos Expansions and Local Sensitivities[END_REF][START_REF] Eldred | Recent Advances in Non-Intrusive Polynomial Chaos and Stochastic Collocation Methods for Uncertainty Analysis and Design[END_REF].

Several UQ methods have been developed with the objective of reducing the number of solutions required to obtain a statistical characterization of the quantity of interest, such as Sparse Grid techniques [START_REF] Smolyak | Quadrature and interpolation formulas for tensor products of certain classes of functions[END_REF] or adaptive mesh generation [START_REF] Foo | The multi-element probabilistic collocation method (ME-PCM): Error analysis and applications[END_REF][START_REF] Witteveen | An adaptive Stochastic Finite Elements approach based on NewtonCotes quadrature in simplex elements[END_REF][START_REF] Witteveen | Refinement Criteria for Simplex Stochastic Collocation with Local Extremum Diminishing Robustness[END_REF][START_REF] Abgrall | A one-time truncate and encode multiresolution stochastic framework[END_REF][START_REF] Tryoen | Adaptive anisotropic spectral stochastic methods for uncertain scalar conservation laws[END_REF]. These techniques can lead to a dramatical reduction of the quadrature points for moderate dimensional problem, provided that the function has some regularity properties. Classical sparse grids [START_REF] Smolyak | Quadrature and interpolation formulas for tensor products of certain classes of functions[END_REF] are constructed from tensor products of one-dimensional quadrature formulas. Some Galerkin-based methods deal with multi-resolution wavelet expansions [START_REF] Maître | Uncertainty propagation using WienerHaar expansions[END_REF][START_REF] Le Maître | Multi-resolution analysis of Wiener-type uncertainty propagation schemes[END_REF], domain decomposition in the random space [START_REF] Babuška | Galerkin finite element approximations of stochastic elliptic differential equations[END_REF], adaptive h-refinement [START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF] for dealing with arbitrary probability distributions.

Among the collocation-based stochastic spectral methods, in [START_REF] Keese | Numerical methods and Smolyak quadrature for nonlinear partial differential equations[END_REF] the authors proposed the use of sparse grid quadrature for stochastic collocation. Older studies show the errors and efficiency of sparse grid integration and interpolation [START_REF] Novak | Simple cubature formulas with high polynomial exactness[END_REF][START_REF] Barthelmann | High dimensional polynomial interpolation on sparse grids[END_REF], Smolyak constructions based on one-dimensional nested Clenshaw-Curtis rules [START_REF] Novak | Simple cubature formulas with high polynomial exactness[END_REF][START_REF] Novak | High dimensional integration of smooth functions over cubes[END_REF] and the integration error of sparse grids based on one-dimensional Kronrod-Patterson rules [START_REF] Petras | On the Smolyak cubature error for analytic functions[END_REF].

An alternative solution is based on approaches attempting to identify the relative importance of the input uncertainties on the output. If some dimensions could be identified as negligible, they could be discarded in a reduced stochastic problem and better statistics estimations could be achieved with a lower computational cost.

Identifying the most influent parameters requires to determine the uncertain inputs which have the largest impact on the variability of the model output. In literature, Global sensitivity analysis (GSA) aims at quantifying how uncertainty in the input parameters of a model contributes to the uncertainty in its output (see for example [START_REF] Borgonovo | Comparison of global sensitivity analysis techniques and importance measures in PSA[END_REF]) where global sensitivity analysis techniques are applied to probabilistic safety assessment models). GSA classifies the inputs according to their importance on the output variations and it gives a hierarchy of the most important ones.

Traditionally, GSA is performed using methods based on the decomposition of the output variance [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF], i.e. ANOVA. The ANOVA approach involves splitting a multi-dimensional function into its contributions from different groups of subdimensions. In 2001, Sobol used this formulation to define global sensitivity indices [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF], displaying the relative variance contributions of different ANOVA terms. In [START_REF] Rabitz | Efficient input-output model representations[END_REF], the authors introduced two High-Dimensional Model Reduction (HDMR) techniques to capture input-output relationships of physical systems with many input variables. These techniques are also based on ANOVA decompositions.

Several techniques have been developed to compute sensitivity indices at low computational cost [START_REF] Saltelli | Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index[END_REF]. In [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF][START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF][START_REF] Blatman | A comparison of three metamodel-based methods for global sensitivity analysis: GP modelling[END_REF], generalized Polynomial Chaos Expansions (gPC) are used to build surrogate models for computing the Sobol's indices analytically as a post-processing of the PC coefficients. In [START_REF] Foo | Multi-element probabilistic collocation method in high dimensions[END_REF], the authors combine multi-element polynomial chaos with an ANOVA functional decomposition to enhance the convergence rate of polynomial chaos in high dimensions and in problems with low stochastic regularity. In [START_REF] Yanga | Adaptive ANOVA Decomposition of Stochastic Incompressible and Compressible Flows[END_REF], the use of adaptive ANOVA decomposition is investigated as an effective dimension-reduction technique in modeling incompressible and compressible flows with high-dimension of random space. In [START_REF] Blatman | Efficient computation of global sensitivity indices using sparse polynomial chaos expansions[END_REF], sparse Polynomial Chaos (PC) expansions are introduced in order to compute sensitivity indices: a PC-based metamodel (a computationally inexpensive polynomial approximation of the relation between inputs and output) which contains the significant terms whereas the PC coefficients are computed by least-square regression.

Other approaches are developed if the assumption of independence of the input parameters is not valid. For instance, new indices have been proposed in [START_REF] Xu | Extending a global sensitivity analysis technique to models with correlated parameters[END_REF][START_REF] Xu | Uncertainty and sensitivity analysis for models with correlated parameters[END_REF] when a linear correlation exists. In [START_REF] Borgonovo | A new uncertainty importance measure[END_REF], the authors introduce a global sensitivity indicator which looks at the influence of input uncertainty on the entire probability distribution without reference to a specific moment of the output (moment independence) and which can be defined also in the presence of correlations among the parameters. In [START_REF] Caniou | Distribution-based global sensitivity analysis in case of correlated input parameters using polynomial chaos expansions[END_REF], a gPC methodology to address global sensitivity analysis for this kind of problems is introduced, while in [START_REF] Borgonovo | Model emulation and moment-independent sensitivity analysis: An application to environmental modelling[END_REF], a numerical procedure is proposed for moment-independent sensitivity methods.

The ANOVA-based analysis create a hierarchy of most important input parameters for a given output when variance is chosen as metrics. A strong limitation of this approach is the fact that the variance can not be considered a general indicator for a complete description of output variations. Roughly speaking, the variance measures the spread of a set of realizations for a random variable, therefore the information related to the direction of this spreading are completely missed. As a result, computation of Higher-Order (HO) statistics is of primary importance, for example the normalized third order, the skewness measures of the non-symmetry of the distribution or the normalized fourth order, the kurtosis measures whether the distribution is peaked or flat. Depending on the problem, a n-order statistics breakdown or, in the following decomposition, could be of interest. Moreover, it seems of primary importance to collect the set of interactions obtained from n-order statistics decomposition for a correct ranking of all the contributions.

For computing HO statistics, the most diffused methods are related to Monte Carlo and quasi-Monte Carlo approaches. Very few papers exist showing the application of polynomial-chaos techniques to the computation of HO statistics [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF][START_REF] Sudret | Uncertainty propagation and sensitivity analysis in mechanical models Contributions to structural reliability and stochastic spectral methods[END_REF].

First objective of this paper is to provide a general method in order to compute the decomposition of high-order statistics, then to formulate an approach similar to ANOVA but for skewness and kurtosis. One of the main results of the paper is to show how skewness and kurtosis can be written mimicking the variance decomposition in term of conditional statistics, i.e. an additive functional relationship always holds. The idea is to compute the most influential interactions not only for the variance but also for higher statistics permitting to improve the sensitivity analysis. This is a fundamental step in order to formulate also innovative optimization methods for obtaining very robust designs by taking into account a complete description of the output statistics. Second objective is to illustrate the correlation between the high-order functional decomposition and the PC-based techniques, thus displaying how to compute each term from a numerical standpoint. Moreover, two reduction strategies, for the resulting polynomial metamodel, are considered. A classical truncation which neglects some non-significant stochastic dimensions and a reduction based on the analysis of the order of interactions. These strategies are evaluated with respect to the complete non-reduced model in terms of their effect on the probability density function of the output.

Several numerical test cases are proposed to demonstrate the interest of the proposed approach. Among them a thermodynamically complex flow in a ORCs turbine cascade, characterized by a significant uncertainty on the physical parameters and on the operating conditions at the turbine inlet [START_REF] Congedo | Numerical investigation of dense-gas effects in turbomachinery[END_REF], is analyzed. Three sources of uncertainties [START_REF] Congedo | Shape optimization of an airfoil in a BZT flow with multiple-source uncertainties[END_REF] are taken into account, namely the thermophysical properties, the inlet boundary conditions and the geometrical parameters of the blade.

The remaining paper is organized as follows. In Section 2, functional decomposition for variance, skewness and kurtosis are presented. In Section 3, the correlation between the functional decomposition and a Polynomial Chaos framework is depicted. Section 4 extends the sensitivity index definitions to high-order statistics and the reduction strategies for a polynomial metamodel are also presented. Several results, showing how the Polynomial Chaos expansion can be used practically to compute high-order statistics and the importance of considering skewness and kurtosis sensitivity indices when ranking the uncertainties/interactions, are presented in Section 5. Moreover in Section 5, the sensitivity analysis of a ORCs turbine is also presented. This numerical investigation aims to display the interest of the proposed approach on a real engineering problem. Eventually, conclusions and future perspectives are drawn in Section 6.

Functional decomposition

Let us consider a real function

f = f ( #» ξ ) with #» ξ a vector of independent and identically distributed random inputs #» ξ ∈ Ξ d = Ξ 1 × • • • × Ξ n (Ξ d ⊂ R d ) and #» ξ ∈ Ξ d -→ f ( #» ξ ) ∈ L 4 (Ξ d , dP).
In the following, we will assume the existence of a density function for the probability measure dP = p(

#» ξ )d #» ξ , where p( #» ξ ) = d i=1 p(ξ i
) is the joint probability density function. We note here that the Sobol' ANOVA decomposition, as well as the Polynomial Chaos approach, requires f ( #» ξ ) ∈ L 2 (Ξ d , dP), however we are interested in defining statistics at least up to the fourth order. Moreover, the random inputs are assumed i.i.d. only for simplicity of exposure.

Let us introduce the definition of the Sobol functional decomposition [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF] for the function f as

f ( #» ξ ) = N i=0 f #» m i ( #» ξ • #» m i ), (1)
where the generic multi-index #» m, of cardinality card( #» m) = d, can contain only elements equal to 0 or 1. The total number of admissible multi-indices #» m i is N + 1 = 2 d ; this number represents the total number of contributions up to the dth-order of the stochastic variables #» ξ . The scalar product between the stochastic vector #» ξ and #» m i is a shorthand employed to identify the functional dependency of f #» m i . In the following, the multi-index #» m 0 = (0, . . . , 0), is associated to the mean term f 0 = E( f ). The remaining N multi-indices are ordered as follows

#» m 1 = (1, 0, . . . , 0) #» m 2 = (0, 1, . . . , 0) . . . #» m d = (0, . . . , 1) #» m d+1 = (1, 1, 0, . . . , 0) #» m d+2 = (1, 0, 1, 0, . . ., 0) . . . #» m N = (1, . . . , 1).
(

) 2
This ordering is assumed for convenience and does not affect in any way the successive ANOVA functional decomposition.

The decomposition (1) is of ANOVA-type in the sense of Sobol [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF] if all the members in Eq. ( 1) are orthogonal with respect to the joint pdf

Ξ d f #» mi ( #» ξ • #» m i ) f #» m j ( #» ξ • #» m j )dP = 0 with #» m i #» m j , (3)
and for all the terms f #» m i , except f 0 ( f 0 is used in the following as shorthand for f #» m 0 ), it holds

Ξ j f #» m i ( #» ξ • #» m i )p(ξ j )dξ j = 0 with ξ j ∈ #» ξ • #» m i . (4)
Each term f #» m i of (1) can be expressed as follows

f #» m i ( #» ξ • #» m i ) = Ξ d-nnz i f ( #» ξ )p( #» ξ ′ i )d #» ξ ′ i - #» m j ⊂ #» m i f #» m j ( #» ξ • #» m j ), (5)
where nnz i is used as shorthand to count the number of non-null elements in #» m i and #» ξ

′ i represents the complement of #» ξ • #» m i to #» ξ , i.e. ( #» ξ • #» m i ) ∪ #» ξ ′ i = #» ξ .
Note that #» m 0 has null cardinality, therefore is always included in the series at the right hand side of Eq. ( 5).

Hereinafter, in order to substantially reduce the complexity of the notation, the integrals are written with respect the probability measure (relative to the multi-index #» m i ):

dP i = p( #» ξ • #» m i )d( #» ξ • #» m i ), (6)
where multiple, two or three, subscripts are adopted to denote the joint probability measure obtained by tensorization of the measures relative to each multi-indices, as for instance dP i j = dP i dP j or dP i jk = dP i dP j dP k . The functional decomposition (1) is usually employed [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF] to compute the contribution of each term to the overall variance, as shown in the next section.

Variance decomposition

ANOVA analysis is based on the variance decomposition in its conditional contributions. Variance can be written in terms of conditional expectation of f and f 2 as

σ 2 = E( f 2 ) -E( f ) 2 . ( 7
)
Starting from Eq. ( 1), it is easy to compute

f 2 ( #» ξ ) = N i=0 f 2 #» m i ( #» ξ • #» m i ) + 2 N i=0 N j=i+1 f #» m i ( #» ξ • #» m i ) f #» m j ( #» ξ • #» m j ). (8)
If the equation ( 8) is integrated in the stochastic space and the orthogonality property (3) is applied, variance can be decomposed as

σ 2 = N i=1 Ξ d f 2 #» m i ( #» ξ • #» m i )dP i = N i=1 Ξi f 2 #» m i ( #» ξ • #» m i )dP i , (9)
where the symbol Ξi is employed to indicate Ξ nnzi for brevity and

E( f ) = f 0 .
Hence, variance can be expressed as the summation of all the conditional contributions

σ 2 = N i=1 σ 2 #» m i , (10)
where

σ 2 #» m i = Ξi f 2 #» m i ( #» ξ • #» m i )dP i . (11)
In the following, a similar decomposition is obtained for the skewness and kurtosis.

Skewness decomposition in conditional terms

In the case of skewness, the major drawbacks is the presence of an higher number of terms to compute with respect to the variance case, i.e. non-null mixed terms exist which, in the case of the variance, are zero due to orthogonality. The first step in order to obtain the skewness decomposition of the function f ( #» ξ ) is to raise the ANOVA series to the third power by employing the multinomial theorem and integrating over the space Ξ. Reducing the number of terms, exploiting the orthogonality of the basis, is possible and a more compact final expression can be obtained (the full derivation is reported in Appendix A)

µ 3 = N p=1 Ξp f 3 #» m p dP p + 3 #» m p #» m q ⊂ #» m p Ξpq f 2 #» m p f #» m q dP pq + 6 N p=1 N q=p+1 N r=q+1 #» m pq = #» m r Ξpq f #» m p f #» m q f #» m r dP pqr . ( 12
)
where Ξi j = Ξ nnz i j and Ξi jk = Ξ nnz i jk . The notation is also simplified by omitting the explicit dependence of the function f #» m i with respect to its coordinates, i.e.

f #» m i = f #» m i ( #» ξ • #» m i ).
Here, a special notation is introduced in order to compute multi-indices as #»

m ab•••z as follows #» m ab•••z = #» m a ⊞ #» m b ⊞ • • • ⊞ #» m z =         m a 1 + m b 1 + • • • m z 1 max 1, z i=a m i 1 , . . . , m a d + m b d + • • • + m z d max 1, z i=a m i d         . (13)
A useful outcome of this decomposition is the possibility to identify the conditional terms related to each single variable or group of variables as expressed for the variance by means of relation [START_REF] Witteveen | An adaptive Stochastic Finite Elements approach based on NewtonCotes quadrature in simplex elements[END_REF]. In the case of skewness, the conditional terms have a more complex expression (except the first order terms, i.e. the terms related to the single variables). This complexity arises from the presence of mixed contributions. To obtain the additive form (as it holds for the variance)

µ 3 = N i=1 µ 3 #» m i , (14)
it is mandatory to identify all the set of indices whose interactions belong to an assigned multi-index #» m i . By considering that to each multi-index #» m i is associated a set of 2 | #» m i | -1 sub-interactions and denoting this set as P i (for instance if #» m i = (1, 1) then the set P i = {(1, 0), (0, 1), (1, 1)}), from the equation [START_REF] Witteveen | Refinement Criteria for Simplex Stochastic Collocation with Local Extremum Diminishing Robustness[END_REF] it is possible to identify each contribution as follows

µ 3 #» m i = Ξi f 3 #» m i dP i + 3 Ξi f 2 #» m i #» m q ∈P i, f #» m q dP i + 6 #» m p ∈P i, #» m p #» m q ∈P i, #» m pq = #» m i Ξi f #» m i f #» m p f #» m q dP i , (15)
where the following shorthand is used

P i, = P i \ { #» m i }.
It is easy to show that Eq. ( 15) reduces to a sum of zero contributions in the case of a function f ( #» ξ ) with Gaussian distribution; for instance assuming f ( #» ξ ) = 3 i=1 ξ i where ξ i ∼ N(0, 1/3), from the convolution of the pdfs we know f ( #» ξ ) ∼ N(0, 1), and µ 3 = 0 with µ 3 #» m i = 0 ∀ #» m i . Note that Eq. ( 15) is also explicitly obtained in the Appendix A.

Kurtosis decomposition in conditional term

The decomposition of the kurtosis is presented in the following. The decomposition based on the functional form Eq. ( 1), after the application of the multinomial theorem and integration, leads to a final expression (for more details see the Appendix B) which reads

µ 4 = N p=1 Ξp f 4 #» m p dP p + 4 #» m p #» m q ⊂ #» m p Ξpq f 3 #» m p f #» m q dP pq + 6 N p=1 N q=p+1 Ξpq f 2 #» m p f 2 #» m q dP pq + 12 N p=1 N q=1 q p N r=q+1 #» m qr \∩ qr ⊆ #» m p Ξpqr f 2 #» m p f #» m q f #» m r dP pqr + 24 N p=1 N q=p+1 N r=q+1 N t=r+1 #» m pq \∩ pq ⊆ #» m rt #» m rt ⊆ #» m pq ⊞∩ rt Ξpqrt f #» m p f #» m q f #» m r f #» m t dP pqrt . (16)
Note that the operator of subtraction by set is employed with the standard notation \, while ∩ rt is a shorthand which indicates the intersection between the two multi-indices #» m r and #» m t , i.e. it contains the variables which belong to both the multi-indices.

As for the skewness, the additive form for the kurtosis

µ 4 = N i=1 µ 4 #» m i , (17)
can be obtained when each specific multi-index is provided. Indeed, the conditional expression for µ 4 #» m i is equal to (see the Appendix B for more details)

µ 4 #» m i = Ξi f 4 #» m i dP i + 4 Ξi f 3 #» m i #» m q ∈P i, f #» m q dP i + 6 #» m p ∈P i #» m p #» m q ∈P i #» m pq = #» m i Ξi f 2 #» m p f 2 #» m q dP i + 12 #» m p #» m p #» m q ∈P i #» m r ∈P i ,r>q #» m p ⊞∩ qr = #» m i Ξi f 2 #» m p f #» m q f #» m r dP i + 24 #» m p ∈P i #» m q ∈P i ,q>p #» m r ∈P i ,r>q t>r, #» m r ∈P i #» m i ⊆ #» m pq ⊞∩ rt #» m i ⊆ #» m rt ⊞∩ pq Ξi f #» m p f #» m q f #» m r f #» m t dP i . (18)

Correlation with Polynomial Chaos Framework

This section is devoted to illustrate how variance, skewness and kurtosis functional decompositions can be computed within the polynomial chaos framework. In this context, an approximation f of the function f is provided

f ( #» ξ ) ≈ f ( #» ξ ) = P k=0 β k Ψ k ( #» ξ ), ( 19
)
where P is computed according to the order of the polynomial expansion n 0 and the stochastic dimension of the problem d

P + 1 = (n 0 + d)! n 0 !d! . ( 20
)
Each polynomial Ψ k ( #» ξ ) of total degree n o is a multivariate polynomial form which involves tensorization of 1D polynomials by using a multi-index

#» α k ∈ N d , with d i=1 #» α k i ≤ n 0 : Ψ k ( #» ξ • #» m ⋆,k ) = d i=1 ψ α k i (ξ i ), (21)
where the multi-index

#» m ⋆,k = #» m ⋆,k ( #» α k ) ∈ N d is a function of #» α k : #» m ⋆,k = (m ⋆,k 1 , . . . , m ⋆,k d ), with m ⋆,k i = α k i / max 1, α k i .
It is important to note here that the number of multi-indices #» α k is P + 1, see Eq. [START_REF] Barthelmann | High dimensional polynomial interpolation on sparse grids[END_REF], which is function of the polynomial degree n 0 . It follows that two multi-indices #» m ⋆,i ( #» α i ) and #» m ⋆, j ( #» α j ) can be equal even if #» α i #» α j , i.e. they contain the same variables but the corresponding 1D polynomials are raised to a different power. By contrast, the ANOVA multi-indices are 2 d and the following condition always holds: #» m i #» m j for i j. For each polynomial basis, ψ 0 (ξ i ) = 1 and then Ψ 0 ( #» ξ ) = 1. Hence, the first coefficient β 0 is equal to the expected value of the function, i.e. E( f ). The polynomial basis is chosen according to the Wiener-Askey scheme in order to select orthogonal polynomial terms with respect to the probability density function p( #» ξ ) of the inputs. Thanks to the orthogonality, the following relation holds

Ξ d Ψ i ( #» ξ )Ψ k ( #» ξ )dP = δ ik Ψ i ( #» ξ ), Ψ i ( #» ξ ) , (22)
where •, • indicates the inner product and δ i j is the Kronecker delta function.

The orthogonality can be used to compute the coefficients of the expansion in a non-intrusive PC framework as follows

β k = f ( #» ξ ), Ψ k ( #» ξ ) Ψ k ( #» ξ )Ψ k ( #» ξ ) ∀k, (23)
where the integrals can be computed by means of any quadrature formula as described in [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF].

Variance decomposition

The second order central moment E( f 2 ) can be evaluated as

E( f 2 ) = Ξ d f ( #» ξ ) 2 dP = Ξ d        P k=0 β k Ψ k ( #» ξ )        2 dP. ( 24
)
This term can be computed exploiting the orthogonality

E( f 2 ) = Ξ d        P k=0 β k Ψ k ( #» ξ )        2 dP = P k=0 β 2 k Ψ 2 k ( #» ξ ) . (25)
As a consequence, variance is

σ 2 = E( f 2 ) -E( f ) 2 = P k=1 β 2 k Ψ 2 k ( #» ξ ) . (26)
Finally, an explicit connection between the last expression and the Eq. ( 9) is found. Thus, each conditional term can be computed as

σ 2 #» m i = k∈K #» m i β 2 k Ψ 2 k ( #» ξ ) , (27)
where K #» m i represents the set of indices associated to the variables included in the vector ( #» ξ • #» m i ):

K #» mi = k ∈ {1, . . . , P} | #» m ⋆,k = #» m ⋆,k ( #» α k ) = #» m i (28)

Skewness decomposition

The PC expansion, Eq. ( 19), raised to the third power by applying the multinomial theorem, and integrating over the whole stochastic space Ξ leads to a form which includes null terms thanks to the orthogonality properties of the PC basis. A compact form (the full derivation is reported in Appendix C) can be written as

µ 3 = P p=1 β 3 p Ψ 3 p ( #» ξ ) + 3 P p=1 β 2 p P q=1 q p β q Ψ 2 p ( #» ξ ), Ψ q ( #» ξ ) ∆ p q + 6 P p=1 P q=p+1 P r=q+1 β p β q β r Ψ p ( #» ξ ), Ψ q ( #» ξ )Ψ r ( #» ξ ) ∆ pqr , (29)
where two functions are introduced for the selection. The first one ∆ p q is defined as follows

∆ p q =        0 if α p j = 0 and m ⋆,q j = 1 1 otherwise ( 30
)
while the function ∆ pqr is defined as

∆ pqr =        0 if m ⋆,p j + m ⋆,q j + m ⋆,r j = 1, 2 1 otherwise. ( 31
)
For a specific #» m i , the previous expression reduces to

s #» m i = p∈K #» m i β 3 p Ψ 3 p ( #» ξ ) +3 p∈K #» m p β 2 p q∈K #» m q #» m pq = #» m i β q Ψ 2 p ( #» ξ ), Ψ q ( #» ξ ) ∆ p q +6 p∈K #» m p q∈K #» m q q≥p+1 r∈K #» m r #» m pqr = #» m i β p β q β r Ψ p ( #» ξ ), Ψ q ( #» ξ )Ψ r ( #» ξ ) ∆ pqr .
(32) Note that the two functions ∆ p q and ∆ pqr should be computed before computing the integral associated to each term for an efficient implementation. Only if this value is one, then the integral need to be truly computed.

Kurtosis decomposition

After the application of the multinomial theorem on Eq. ( 19) and its integration over Ξ, the expression for the kurtosis is obtained, in compact form as (the full form is derived in Appendix D)

µ 4 = Ξ d f ( #» ξ ) -β 0 4 dP = Ξ d         P p=1 β p Ψ p ( #» ξ )         4 dP = P p=1 β 4 p Ψ 4 p ( #» ξ ) + 4 P p=1 β 3 p P q=1 q p β q Ψ 3 p , Ψ q ∆ p q + 6 P p=1 β 2 p P q=p+1 β 2 q Ψ 2 p , Ψ 2 q + 12 P p=1 β 2 p P q=1 q p β q P r=q+1 r p β r Ψ 2 p , Ψ q Ψ r ∆ p qr + 24 P p=1 P q=p+1 P r=q+1 P t=r+1 β p β q β r β t Ψ p Ψ q , Ψ r Ψ t ∆ pqrt , (33)
where the function ∆ p q is already introduced in (30), while the others two functions are defined as follows

∆ p qr =        0 if α p j = 0 and m ⋆,q j + m ⋆,r j = 1, 2 1 otherwise and ∆ pqrt =        0 if m ⋆,p j + m ⋆,q j + m ⋆,r j + m ⋆,t j = 1, 2 1 otherwise. ( 34
)
The functions ∆ q p , ∆ pqr , ∆ p qr and ∆ pqrt can be computed efficiently before evaluating the integral terms of the expansion. A detail of a possible algorithm for their evaluation is reported in Appendix E.

In the case of the conditional contribution associated to a specific multi-index #» m i , it holds

µ 4 #» m i = k∈K #» m i β 4 k Ψ 4 k ( #» ξ ) + 4 p∈K #» m p β 3 p q∈K #» m q -{p} #» m p ⊞ #» m q = #» m i β q Ψ 3 p , Ψ q ∆ p q + 6 p∈K #» m p β 2 p q∈K #» m q -{p} #» m p ⊞ #» m q = #» m i β 2 q Ψ 2 p , Ψ 2 q +12 p∈K #» m p β 2 p q∈K #» m q -{p} β q r∈K #» m r r≥q+1 #» m pqr = #» m i β r Ψ 2 p , Ψ q Ψ r ∆ p qr + 24 p∈K #» m p q∈K #» m q q≥p+1 r∈K #» m r r≥q+1 t∈K #» m t t≥r+1 #» m pqrt = #» m i β p β q β r β t Ψ p Ψ q , Ψ r Ψ t ∆ pqrt . ( 35
)

Sensitivity indices and truncation error estimation

As introduced by Sobol [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF], sensitivity indices for variance can be computed for each conditional contribution following Eq. ( 10):

σ 2,SI #» m i = σ 2 #» m i σ 2 . ( 36
)
Here, we introduce additional sensitivity indices based on the decomposition of skewness and kurtosis 1 . From the definition of the conditional terms in [START_REF] Maître | Uncertainty propagation using WienerHaar expansions[END_REF] and (18) it follows

s SI #» m i = s #» m i s k SI #» m i = k #» mi k . (37)
In the case of the total sensitivity index (TSI) it is necessary to compute the overall influence of a variable. This sensitivity index can be computed summing up all the contributions in which the chosen variable is present

TSI j = ξ j ∈( #» ξ • #» m i ) σ 2,SI #» m i TSI s j = ξ j ∈( #» ξ • #» m i ) s SI #» m i TSI k j = ξ j ∈( #» ξ • #» m i ) k SI #» m i . (38)
Moreover, following [START_REF] Caflisch | Valuation of Mortgage Backed Securities Using Brownian Bridges to Reduce Effective Dimension[END_REF], it is possible to define two distinct effective dimensions for the ANOVA representation (1), namely the effective dimension δ tr in the truncation sense

f ( #» ξ ) = #» m i f #» m i ≃ 0<i≤δ tr f #» m i = f tr , (39)
and the effective dimension δ sup in the superposition sense

f ( #» ξ ) = #» m i f #» m i ≃ nnz i ≤δ sup f #» m i = f sup . (40)
The two integers δ t and δ s represent the number of terms which are needed in order to reach the required amount of variance, i.e. σ2 ( f tr ) ≤ C tr σ 2 ( f ) and σ 2 ( f sup ) ≤ C sup σ 2 ( f ) respectively. Generally, constants C tr and C sup are very close to the unity; in [START_REF] Caflisch | Valuation of Mortgage Backed Securities Using Brownian Bridges to Reduce Effective Dimension[END_REF] both are chosen to be 0.99. It is important to note here that the dimension δ tr is dependent on the ordering of the terms employed in [START_REF] Xiu | The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations[END_REF]. From a practical point-of-view it is common to approximate the function of interest with a model including only first order interaction terms, i.e. terms which depend only on one single variable 2 . Hereinafter, the generic term metamodel is used to denote an approximation of the true function; it can assume any functional form when not differently specified.. This case is equivalent to choose

δ sup = 1, thus nnz i ≤1 f #» mi = f I .
When the superposition dimension is unitary, i.e. only the terms depending on each single variable are retained, it is possible to deduce explicitly the central moments corresponding to this metamodel

Ξ d ( f ( #» ξ ) -f 0 ) n dP ≈ Ξ d ( f I ( #» ξ ) -f 0 ) n dP = Ξ         nnz i =1 f #» m i         n dP. ( 41
)
For n = 2 (variance) from the orthogonality of the ANOVA terms it follows

Ξ d ( f I ( #» ξ ) -f 0 ) 2 dP = nnz i =1 Ξi f 2 #» m i dP i = σ 2 I = C σ 2 σ 2 , with C σ 2 ≤ 1. ( 42
)
The case of the third order is obtained with n = 3

Ξ d ( f I ( #» ξ ) -f 0 ) 3 dP = nnz i =1 Ξi f 3 #» m i dP i = µ 3 I = C s µ 3 . (43)
In this case the constant C s can be either negative or zero. It is important to note that the terms of the kind Ξi j f 2 #» m i f #» m j dP i j , with i j, are null for orthogonality. The skewness associated to the metamodel including only first order interaction terms has the same functional form of the variance, i.e. it is additive with respect to the contributes relates to each single variable.

The kurtosis results in a more complex form. By raising to the fourth power, n = 4, it follows

Ξ d ( f I ( #» ξ ) -f 0 ) 4 dP = nnz i =1 Ξi f 4 #» m i dP i + 6 nnz i =1 nnz j =1 j>i Ξi j f 2 #» m i f 2 #» m j dP i j = µ 4 I = C k µ 4 , with C k ≥ 0. ( 44
)
The fourth central moment µ 4 I , in this case, includes the contributions associated to the second order interactions. Thus any of the terms

f 2 #» m i f 2 #» m j dP i j belongs to k #» m i j , where #» m i j = #» m i ⊞ #» m j .
The kurtosis cannot be obtained only summing the first order contributions. Hence, this metamodel exhibits an higher kurtosis than the one corresponding only to the sum of its first order interaction terms.

If the usual definitions of skewness or kurtosis are needed, it is necessary to normalize the corresponding central moments by a factor, namely the third and fourth powers of the standard deviation. Hence, the error on the second central moment, i.e. the variance, needs to be taken into account as follows

s I = µ 3 I (σ 2 I ) 3/2 = C s (C σ 2 ) 3/2 µ 3 (σ 2 ) 3/2 k I = µ 4 I (σ 2 I ) 2 = C k (C σ 2 ) 2 µ 4 (σ 2 ) 2 (45)
The estimations obtained in this section allow to examine the quality of a metamodel including only first order interaction terms by using additional information with respect the approach based only on the variance. In Section 5 these estimations are shown to be very useful to avoid situation in which a metamodel is accepted only based on the variance, while in reality it produces poor results in terms of high-order statistics. For instance, it is very common to obtain pdf which cannot preserve the asymmetry (skewness) properties of the original probability density function. Moreover, the previous estimations [START_REF] Guardone | Maximum intensity of rarefaction shock waves for dense gases[END_REF] are only based on the computation of the first order conditional terms (see [START_REF] Maître | Uncertainty propagation using WienerHaar expansions[END_REF] and [START_REF] Caniou | Distribution-based global sensitivity analysis in case of correlated input parameters using polynomial chaos expansions[END_REF] with nnz i = 1), thus their evaluation results to be always affordable numerically.

Moreover, the previous estimations can be adopted to select the best possible metamodel, in term of polynomial degree given the simulations at disposal. A classical non-intrusive implementation of a PC method starts from the knowledge of n = 1, . . . , N functional evaluations of the model. For each of these realization a weight w n is associated according to the quadrature rule desired (it can be either a full tensorization of Gauss-based 1D quadrature rules, Sparse Grids or other more sophisticated choices, see also [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF]). The pseudocode which follows is an example of a possible use of the estimations reported in Eq. 45 for the selection of the optimal polynomial degree of a metamodel including only first order interaction terms: 3 , and µ 4

n 0 = 1, iter = 1 and ε = user defined tolerance while max(∆C σ 2 , ∆C s /(C σ 2 ) 3/2 , ∆C k /(C σ 2 ) 2 ) > ε & n 0 > 1 do for i = 1, . . . , d do for k ∈ K #» m i do β k = f ( #» ξ ),Ψ k ( #» ξ ) Ψ k ( #» ξ )Ψ k ( #» ξ ) end end for n = 1, . . . , N do f (n) #» m i = k∈K #» m i β k Ψ (n) k end for i = 1, . . . , d do σ 2 I,i = N n=1 w n ( f (n) #» m i ) 2 , µ 3 I,i = N n=1 w n ( f (n) #» m i )
I,i = N n=1 w n ( f (n) #» m i ) 4 ( 46
)
for j = i + 1, . . . , d do µ 4 I,i = µ 4 I,i + 6 N n=1 w n ( f (n) #» m i ) 2 ( f (n) #» m j ) 2 end end
Assembling the central moments:

σ 2 I = d i=1 σ 2 I,i , µ 3 I = d i=1 µ 3 I,i and µ 4 I = d i=1 µ 4 I,i (47)
Evaluation of the central moments by collocation: 3 , and

µ = N n=1 w n f (n) , σ 2 = N n=1 w n ( f (n) ) 2 -(µ) 2 , µ 3 = N n=1 w n ( f (n) ) 3 -3µσ 2 -(µ)
µ 4 = N n=1 w n ( f (n) ) 4 -4µ 3 µ -6(µ) 2 σ 2 -(µ) 4 . (48)
Evaluation of the constants for the first order interaction metamodel:

C σ 2 = σ 2 I σ 2 , C s = µ 3 I µ 3 and C k = µ 4 I µ 4 (49) if iter > 1 then ∆C σ 2 = C σ 2 | iter -C σ 2 | iter-1 ∆C s = C s (C σ 2 ) 3/2 iter - C s (C σ 2 ) 3/2 iter-1 ∆C k = C k (C σ 2 ) 2 iter - C k (C σ 2 ) 2 iter-1 (50) end n 0 = n 0 + 1 and iter = iter + 1 end
Algorithm 1: Pseudocode for the choice of the optimal order for the first order interactions metamodel.

This algorithm requires only n 0 d coefficients out of the P+1 = (n 0 +d)!/(n 0 !d!) needed for a polynomial expansion of total order n 0 . Indeed, the set K #» m i for a multi-index #» m i with nnz i = 1 reduces to a set of n 0 elements, i.e. one for each polynomial order up to the total order n 0 . Therefore, without requiring additional model evaluations the total degree n 0 can be raised up to reaching the convergence of the constants C σ 2 , C s and C k . The obtained polynomial degree represents the optimal degree for a metamodel including only first order interactions. In such a case, the metamodel obtained is the one including the minimum number of coefficients, however it is important to note that in complex models the metamodel including only first order contributions can furnish misleading interpretations of the tail of the distributions. In a case like that, the knowledge of the high-order moments decomposition is sufficient to evaluate, a priori, the quality of the metamodel. Relevant examples for this situation are reported in Section 5.5.

Numerical results

In this section, the importance of considering metrics based on high-order statistics for global sensitivity analysis is demonstrated through some numerical examples and a complex fluid flow application. The numerical test section is organized as follows. In Section 5.1, the previous relations for the computation of the high-order conditional terms are demonstrated numerically by showing the convergence properties of PC with respect to the analytical conditional high-order statistics. In Section 5.2, a comparison is performed between the information obtained by an analysis based on the variance or on high-order conditional contributions for numerical test problems with different kind of interactions between parameters. How reducing the model dimension in the truncation and in the superposition sense is described in Sections 5.3 and 5.4, respectively. Finally, in Section 5.5, the high-order decomposition analysis of a complex flow in a Organic Rankine Cycles (ORCs) turbine is presented, thus highlighting the importance of this kind of study during the design process of a complex system.

Computing conditional statistics by means of PC

The high-order statistics (the values of variance, skewness or kurtosis) can be computed using the same set of deterministic evaluation of the models, i.e. the same number of functional evaluations f = f ( #» ξ ) in the same sample points, by evaluating the first order coefficients of the PC expansion for f , f 2 , f 3 and f 4 corresponding to the expected values of the four functions. Hereinafter, this approach is referred to as collocation. It is important to remark that the collocation approach does not provide any kind of metamodel for the function f , nor the possibility to compute conditional terms. Then, this approach is employed only for a comparison with respect to the PC series and for assessing the convergence of the expansion.

Let consider the following function

f ( #» ξ ) = d i=1 sin(πξ i ) (51)
where each variable ξ i ∼ U(0, 1) with an increasing dimension d up to three. In the following, statistical moments as well as sensitivity indices (relative) errors are systematically computed with respect to the analytical solution.

In Figure 1, convergence of the statistical moments is reported as a function of the number of functional evaluations for the dimension d = 2. A number of simulations equal to N = 120 is needed to reach a relative error of order O(10 -4 ) for the kurtosis. Also the collocation approach, as expected, converges faster, but, as already discussed, is limited to the computation of the full central moments. Now, conditional statistics can be computed using a PC approach using Eqs. ( 29) and [START_REF] Xu | Uncertainty and sensitivity analysis for models with correlated parameters[END_REF]. In Figure 2, we show first-order interaction terms σ 2 1 , µ 3 1 and µ 4 1 (where for symmetry

σ 2 1 = σ 2 2 , µ 3 1 = µ 3 2 and µ 4 1 = µ 4
2 ) and second order interaction terms (σ 2 12 , µ 3 12 , µ 4 12 ) errors computed with respect to the analytical solution. Statistics are well converged at N = 120. Then, the case with d = 3 is considered. In Figures 3 and4, convergences of statistical moments and conditional statistics are reported, respectively. Convergence, for both statistical moments and conditional statistics with respect to their exact analytical counterparts, is attained at nearly N = 1500.

High-order indices analysis for global Sensitivity Analysis

The importance of including high-order conditional terms in the analysis is demonstrated in this section by means of several model functions. Let us consider, the classical Sobol function (four dimensions)

f ( #» ξ ) = 4 i=1 |4ξ i -2| + a i 1 + a i , (52)
where ξ i ∼ U(0, 1). Two possible choices of the coefficients a i are considered

• a i = (i -1)/2, namely linear g-function f glin ;

• a i = i 2 , namely quadratic g-function f gquad .

In Figure 5, Sensitivity Indices (SI) for the linear g-function f glin are reported. Several differences can be observed between the sensitivity indices computed on the variance or on high-order moments. The variance-based ranking illustrates that the first-order sensitivity indices, i.e. the indices relative to the first order interaction terms are higher than the second order ones, while these last ones are higher than the third and fourth order ones. This is not the case for skewness and kurtosis, where the second-order are the larger contributions. This behavior reveals that the variance is able to capture the ranking of the variables, but not the relative importance associated to higher-order interactions between variables. From a practical point of view, this error can lead to wrong decisions in a dimension reduction strategy as it will be shown in Sections 5.3 and 5.4. Quantitatively, the variance based only on first-order contributions exceeds 0.8, while skewness and kurtosis do not attain 0.1. In Table 1, the total sensitivity indices for the four variables are reported. It is evident that the ranking of variables is not influenced by the statistical moment chosen as metrics, but their relative importance is considerably different. The same functional form, Eq. ( 52), can lead to slightly different results if the quadratic function coefficients are considered. In Figure 6, the sensitivity indices for the g-function with a quadratic dependence of the coefficients are reported. In this case, the difference between the first order contributions and high-order terms is even more evident. By considering the variance, first-order contributions exceed 0.98, while a value larger than 0.5 is computed for highorder interactions for skewness or kurtosis. In this case, the contribution of the first variable exceeds 0.8, while for the skewness or kurtosis, in order to attain the same level, it is necessary to include also the contributions related to the interaction between the first and second variable. In Table 2, total sensitivity indices are reported for the four variables. In this case, variance contributions for both the third and fourth variables are below 0.05, while for both skewness and kurtosis, only the fourth variable contribution takes a TSI value of 0.04, which could be considered low.

Variable

A low level of TSI for the variables ξ 3 and ξ 4 could suggest to reduce the model to the first two variables or neglect the contributions related to the order higher than one in a polynomial metamodel, as for instance the truncated PC series. Let consider the following functions:

f 1 = ξ 1 exp       ξ 2 ξ 2 3 + 1       + ξ 1 ξ 2 f 2 = 3 i=1 2ξ i + 1 2 , (53)
where the parameters are ξ i ∼ U(0, 1). Sensitivity indices associated to the first function f 1 are reported in Figure 7. For function f 1 , the most important variable is ξ 1 . For the variance, the first-order sensitivity index relative to ξ 1 is also the most important SI. On the contrary, for both skewness and kurtosis, the highest SI is associated to the second-order interaction between the first and the second variable. In this case, the inspection of the total sensitivity indices, reported in Table 3, suggests that the third variable ξ 3 is meaningless with respect to the variance. The TSI associated to ξ 3 are lower than the limit proposed in [START_REF] Gao | Efficient solution of ordinary differential equations with high-dimensional parametrized uncertainty[END_REF] to identify a negligible uncertainty. However, if this information is used together with the high-order total sensitivity indices information, the choice of neglecting the third variable should be considered more carefully. The results of a model reduction decision, totally based on variance measures, is further discussed in the following section.

Function f 2 is reported here to underline the difference between the measure of sensitivity associated to the variance and to the higher-order moments. In particular, the functional form of f 2 , Eq. (53), includes an equal contribution of three variables. However, looking at Figure 8, it is possible to note that the variance is concentrated only on firstorder contributions of the single variables and their sum exceeds 0.9. The skewness and kurtosis contributions, on the contrary, are concentrated on second-and third-order interactions, respectively. Even if the sum of the first-order variance contributions exceeds 0.9, a reduction of the model in the superposition sense (i.e. by neglecting the high orders of interaction), could lead to wrong conclusions, as explained in Section 5.4. The skewness associated to a model including only first-order contributions does not include the skewness information about the probability distribution of the output.

Values for the total sensitivity indices are reported in Table 4 for this case. It is interesting to note that the sum of the total sensitivity indices over the three variables is higher for skewness and kurtosis than for the variance. Then, they reveal an intrinsically high-order interactions (see Eq. (53) for f 2 definition). Numerical test cases presented in this section illustrate how information relative to variance-based sensitivity indices are insufficient in order to understand the true dependence of a model from its variables. Moreover, for a function that is known by points, i.e. for example experimental observations or runs of a numerical code, the sensitivity indices on the skewness and on the kurtosis could be very helpful to capture some interactions between subset of variables, much more than the variance.

Variable

This could be even more important if the aim is to build/reduce the dimensionality of the problem and to build an accurate metamodel.

Dimensional reduction in the truncation sense

We start from the quadratic g-function, Eq. ( 52). From the analysis conducted in the previous section (see Table 2), the third and fourth variable seem to be meaningless for the variance-based indices. Their total sensitivity indices sum up to 0.05 for the variance, while exceed 0.15 for both skewness and kurtosis. Considering only the sensitivity indices computed on the variance, one could be tempted to neglect the variables ξ 3 and ξ 4 as follows

f G1 = f 0 + f 1 (ξ 1 ) + f 2 (ξ 2 ) + f 12 (ξ 1 , ξ 2 ) f G2 = f 0 + f 1 (ξ 1 ) + f 2 (ξ 2 ) + f 12 (ξ 1 , ξ 2 ) + f 3 (ξ 3 ) + f 13 (ξ 1 , ξ 3 ) + f 23 (ξ 2 , ξ 3 ) + f 123 (ξ 1 , ξ 2 , ξ 3 ), ( 54
)
where in the first case, f G1 , both variables are neglected; on the contrary for f G2 only ξ 4 is neglected. In this case, the ANOVA terms and the statistics can be computed analytically. In Table 5, the percentage errors, for the first four central moments, are reported with respect to the analytical exact solution for both the reduced models f G1 and f G2 .

Function × 100 errors related to the reduced g-function f G1 and f G2 .

In Table 5, it is evident that an error of only 5% on the variance can correspond to a much larger error on the higher moments. This behavior is confirmed by at Figure 9, where the probability density function is computed for both f G1 and f G2 and compared with the complete function, Eq. (52). In this case, f G1 which includes only the first two variables can not reproduce the tails while a good approximation is attained elsewhere. In this case, the pdf is bounded between 0.4 and 1.8. If the third variable is included in the reduced model, both variance and kurtosis are computed with an error lower than 5%, while the error on the skewness remains lower than 8%. The total sensitivity indices associated to the fourth variable are reported in Table 2 and it is lower than 5% for the three moments. The improvement of the model given by including the third variable, f G2 , is evident in Figure 9, where the pdf of the reduced model better approximates the pdf of the complete function.

From a practical point-of-view, the dimension reduction is commonly accomplished by neglecting some inputs. On the contrary, for an analytical function, it is possible to compute the constant values to choose, for both ξ 3 and ξ 4 , in order to obtain a reduced model that preserves both the expected value and the variance of the original complete model. Of course, both requirements cannot be satisfied at the same time, but a set of values satisfying the mean and the variance can be obtained analytically requiring that

E( f ( #» ξ )) = |4 ξ j -2| + a j 1 + a j = 1 E( f 2 ( #» ξ )) = |4 ξ j -2| + a j 1 + a j 2 = 1 0 |4ξ j -2| + a j 1 + a j 2 dξ j . (55)
The following values can be analytically computed for the two variables: ξ 3 = {1/4, 3/4, 91/120, 29/120} and ξ 4 = {1/4, 3/4, 77/102, 25/102}.

In Figure 10, the pdf associated to the complete quadratic g-function with parameters ξ 3 and ξ 4 neglected is reported. In particular, all the possible choices of ξ 3 and ξ 4 are employed obtaining distinct probability functions which are very close each other. This behavior is an additional proof that the third and fourth variable are not significant and their values can only slightly affect the functional response.

From Figure 10, it is evident that neglecting some parameters in order to assure the correctness of the mean and the variance yields a pdf very close to that one obtained by neglecting entirely the ANOVA terms. From a practical point-of-view, the analysis of the reduced model can be carried out both with the ANOVA reduced model (if it is analytically possible to compute the integrals) or by imposing the parameters satisfying the expected value and variance. However, the main point here is that variance-based sensitivity analysis should be supplemented by highorder sensitivity analysis for building a reduced model which does not deteriorate the pdf of the results especially with respect to the tails.

Dimensional model reduction in the superposition sense

In this section, the problem of the truncation is analyzed from a different perspective, in a so-called superposition sense. This means that the dimension of the model is not reduced in terms of number of variables, but in terms of order of interaction between variables for the polynomial metamodel. Note also that, if the function is approximated by means of a PC series of total degree n 0 , all the interactions of order n 0 + 1 are lost. From a practical standpoint, the error related to the truncation of the PC series is greater than the error of the truncation of the ANOVA approximation at a certain order. This is due to the approximation of each single term of the ANOVA expansion via a truncated polynomial series; the ANOVA functional decomposition contains 2 d terms which are approximated by a finite PC series. In the following, the reduced model is computed analytically. This case relies on a perfect knowledge of the reduced model. However, to obtain an equivalent result for a generic function a non truncated PC series would be necessary.

The first example considered is the linear g-function, Eq. ( 52). Results presented in Section 5.2 in terms of sensitivity indices (see Figure 5) and total sensitivity indices (see Table 1) illustrate some main features: the first order interaction seems to be important enough to entirely approximate the model; the contributions related to the first order interactions exceed 0.8 for the variance, but it is much less for both skewness and kurtosis. Two different reduced models are considered in this case: the first-order model f O1 and the second-order one f O2 , described by the following equations

f O1 = f 0 + f 1 (ξ 1 ) + f 2 (ξ 2 ) + f 3 (ξ 3 ) + f 4 (ξ 4 ) f O2 = f O1 + f 12 (ξ 1 , ξ 2 ) + f 13 (ξ 1 , ξ 3 ) + f 14 (ξ 1 , ξ 4 ) + f 23 (ξ 2 , ξ 3 ) + f 24 (ξ 2 , ξ 4 ) + f 34 (ξ 3 , ξ 4 ). ( 56
)
In Table 6, the contributions relative to each moment, for the two models, are reported, where the models are obtained by a truncated PC series and their exact counterparts, f ex O1 and f ex O2 are computed analytically. It is important to note here that the presence of the absolute value, in both the version of the g-function, prevents the spectral convergence of the PC expansion [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF]. Hence, the difference between the PC and the analytical values in Table 6 depends on the PC series truncation.

Figure 11 illustrates the PDF for the complete and the reduced model. Note that including a large amount of variance could lead to a largely inaccurate metamodel if the information on the variance are not supplemented by those obtained by the analysis of high-order moments. From Table 6, it is evident that even if the variance related to the first order terms exceeds the 80% of the total variance, the corresponding skewness and kurtosis are very low. The situation is evident in the probability density function associated to the reduced model f O1 reported in Figure 11 where the pdf corresponding to f O1 does not contain any skewness (it is perfectly symmetric). The situation greatly improves including contributions up to the second order of interactions between variables. This is supported by the statistics reported in Table 6 where, even if the improvement in terms of variance is reduced, a better approximation of both skewness and kurtosis is achieved.

The second example is the function f 2 , Eq. ( 53). The results reported in Section 5.2 show than the first order terms represent more than 90% of the variance while they correspond to the 0% for the skewness and they contribute to less than 15% for the kurtosis. In this case, looking at the sensitivity indices relative to the variance, a model including only first order interaction terms could appear as a good approximation of the complete function. However, considering the first and the second order defined as follows

f O1 = f 0 + f 1 (ξ 1 ) + f 2 (ξ 2 ) + f 3 (ξ 3 ) f O2 = f O1 + f 12 (ξ 1 , ξ 2 ) + f 13 (ξ 1 , ξ 3 ) + f 23 (ξ 2 , ξ 3 ), ( 57
)
the computation of the pdf reveals the importance of the high order interaction terms. In Figure 12, the pdf for the complete model and the first and second orders are reported. Even if more than 90% of the variance is included in the first order model, its pdf contains no information about the skewness and the tails.

Analysis of a ORCs turbine operating under uncertainty conditions

The last numerical example is a turbine blade of a two dimensional VKI LS-59 cascade, a configuration which has been widely studied [START_REF] Congedo | Numerical investigation of dense-gas effects in turbomachinery[END_REF][START_REF] Kiock | The transonic flow through a plane turbine cascade as measured in four european wind tunnels[END_REF]. An unstructured CFD dense-gas solver is used to ensure the reliability of the computed results for dense gas flows (for more details see Ref. [START_REF] Congedo | Numerical investigation of dense-gas effects in turbomachinery[END_REF]). The two-dimensional flow domain is discretized by a structured C-grid comprised of 192x16 cells (see Figure 13(a)). The boundary conditions are imposed as follows: at the inlet and outlet boundaries, non-reflecting boundaries are applied using the method of characteristics; a slip condition is imposed at the wall, which uses multi-dimensional linear extrapolation from interior points to calculate the wall pressure; periodicity conditions are prescribed at the inter-blade passage boundaries. The siloxane dodecamethylcyclohexasiloxane (C 12 H 36 Si 6 O 6 ), commercially known as D 6 , is the fluid considered in this study. The physical properties of D 6 are reported in Table 7. The Peng-Robinson-Stryjek-Vera model (PRSV) equation is used as thermodynamic model for D6, which is provided in the non-dimensional form (critical point properties are used as reference) as follows:

P r (T r , v r ) = T r /Z c v r -b r - a r v 2 r + 2bv -b 2 , ( 58
)
where a r and b r are substance-specific parameters related to the fluid critical-point properties:

a r = 0.457235/Z 2 c α(T r ), b r = 0.077796/Z c , (59)
and Z c is the critical compressibility factor, i.e. Z c = (P c v c )/(RT c ). It can be computed by enforcing the passage of the critical isotherm through the critical point, and then it is the solution of a cubic equation. The correction factor α is given by

α(T ) = 1 + K(1 -T 1/2 r ) 2 , ( 60
) with K = γ 0 + γ 1 ω -γ 2 ω 2 + γ 3 ω 3 , ( 61
)
where the parameter ω is the fluid acentric factor and the coefficients (γ 0 , γ 1 , γ 2 , γ 3 ) are equal to (0.378893, 1.4897153, -0.17131848, 0.0196554), respectively.

To close the thermodynamics relations, a power law for the ideal gas specific heat at constant volume where n a fluid-dependent parameter and T c the critical temperature. The equation of state for the internal energy, e, is computed by exploiting the compatibility relations (see [START_REF] Cinnella | Sensitivity analysis of dense gas flow simulations to thermodynamic uncertainties[END_REF] for more details). In summary, the PRSV model depends on the fluid acentric factor ω, the ideal-gas specific heat at the critical temperature c v∞ (T c ), and a fluid-dependent parameter n.

(c v,∞ ) is assumed c v,∞ (T ) = c v,∞ (T c ) (T r ) n , (62)
Three sources of uncertainties are considered in this study (globally four uncertainties parameters): uncertainty on the operating conditions, namely the inlet total temperature, T in /T c and inlet total pressure, p in /p c ; uncertainty on the thermodynamic model, namely n (from previous studies, the large predominance of n with respect to ω, c v∞ has emerged [START_REF] Congedo | TSI metamodels-based multi-objective robust optimization[END_REF][START_REF] Congedo | Shape optimization of an airfoil in a BZT flow with multiple-source uncertainties[END_REF]); and uncertainties on geometrical parameters, namely the blade thickness φ. Following previous investigations [START_REF] Colonna | Real-gas effects in organic rankine cycle turbine nozzles[END_REF], 3.0% of uncertainty for the temperature and pressure levels at the inlet conditions is considered. Since the D6 fluid is considered, the parameter n is assumed to be affected by 6% uncertainty (mean values equal to 0.5729) following [START_REF] Cinnella | Sensitivity analysis of dense gas flow simulations to thermodynamic uncertainties[END_REF][START_REF] Colonna | Real-gas effects in organic rankine cycle turbine nozzles[END_REF], with values of ω and c v∞ equal, respectively to 0.7361 and 105.86. An uncertainty of 2% for the thickness φ is considered.

Performance of the turbine cascade can be evaluated by using several output criteria. Statistics decomposition is applied to some specific quantities of interest, which define global performances and efficiency of the turbine. In particular, the following criteria are selected:

• ∆h, which is the enthalpy variation through turbine stage;

• The power output per unit depth (PO) expressed as ∆h • ṁ/w mol [W], where ṁ is the mass flow rate and w mol is the molecular weight;

• the relative temperature variation or Carnot factor ∆T/T inlet , where ∆T represents the variation of temperature between turbine inlet and outlet boundaries;

• the turbine isentropic efficiency ∆h/∆h ideal , where ∆h represents the variation of the static enthalpy between turbine inlet and outlet boundaries, and ∆h ideal is the variation of enthalpy when an ideal isentropic transformation with the same initial conditions and pressure ratio as the real flow.

In Figures 14 and15 the sensitivity indices are reported for the four outputs of interest. The results are based on 4096 simulations obtained at the locations of the tensorization of 1D Legendre polynomials with order 7 in each direction. The PC approximation has been chosen, after a convergence study not reported here for brevity, to use a total polynomial degree of order 5. For all the outputs the third order interaction terms are negligible for any statistic moments. The situation changes when first and second order contributions are considered. The first order terms contribute almost entirely to the overall variance, while for skewness and kurtosis their effect is as important as the second order terms. The interesting feature of the enthalpy variation (Figure 14(b)), the Carnot factor (Figure 15(a)) and the power output (Figure 15(b)) is the presence of negative contributions associated to the skewness. For the enthalpy variation and the power, the first order terms sum up to negative values of the order of the unity. Such a behavior means that the skewness associated to the first order metamodel shows a skewness roughly equal in magnitude to the complete function, but with an opposite sign. The effect of a similar situation is evident comparing the probability density functions of distinct metamodels associated to different superposition lengths.

In Table 8 the a priori linear metamodel estimations (see Section 4) are applied for the turbine case.

In Figures 16 and17 the pdfs for the four outputs of interest are reported. For all the outputs it is evident that the metamodel including up to second order interaction terms is able to represent accurately the original function. For the efficiency, Figure 16(a), and the Carnot factor, Figure 17(a), the skewness associated to the model including only first order interaction terms has the same sign of the overall third moment. For the other outputs, namely the enthalpy variation, Figure 16(b), and the power output, Figure 17(b), the skewness of the first order model has the opposite sign and the effect is an inversion of the tails. However, we note here that the effect of the sign inversion is limited as it is Output associated to small values of skewness for the complete model, namely -8.229 × 10 -2 for the enthalpy variation and -0.1539 for the power output.

µ 2 I C σ 2 µ 3 I C s µ 4 I C k C s /(C σ 2 ) 3/2 C k /(C σ
To further demonstrate the relevance of the high-order statistics, for the right tails of the enthalpy variation and power output, for which the metamodel including only first order interaction terms is providing poor approximations in term of skewness, the probability of exceeding a fixed value is reported in Figure 18. The probability of exceeding a fixed value is obtained as the difference between one and the cumulative density function (CDF) value. For instance, if the probability of exceeding the value of 1.5 for the power output is needed, then the metamodel including first order interactions can predict a probability of order 10 -3 , while the complete model (or the more accurate metamodel) attain a value of order 10 -2 .

In Section 4 a procedure for the optimal selection of the polynomial degree for a metamodel exploiting the information related to the high-order statistics has been introduced. In Figures 19 and20 the convergence of the constants C σ 2 , C s /(C σ 2 ) 3/2 and C k /(C σ 2 ) 2 are reported as function of the number of coefficients of the metamodel, n 0 d, for the different outputs of the ORC turbine.The convergence is reported in term of absolute value of the difference between two successive polynomial degree and it shows how the variance is attained with low errors for low polynomial degrees. In all the cases, but the efficiency, an error of 10 -2 is reached with order n 0 = 2 for the constant relative to the variance, while the convergence is slower for the constants related to the high-order statistics. In particular, if a threshold of 10 -1 is fixed for the algorithm reported in Section 4 it is evident that an higher, n 0 > 2, is always needed. The outputs with slower convergence of the constants related to high-order statistics are the efficiency and the power output, reported in Figures 19(a) and 20(b), respectively. For these outputs, in Figures 21 and22 the probability density function and the cumulative density function are reported. For both efficiency and power output the threshold of ε = 10 -1 is reached with n 0 = 4 for the high-order statistics, while n 0 = 2 is high enough to provide convergence for the constant related to the variance.

Conclusions

We present a decomposition of high-order statistics and the importance of using this information for reducing the complexity of a uncertainty analysis.

We also present a correlation between the functional decomposition, as depicted by Sobol, and the polynomial chaos development which enables to identify each term of the decomposition, drawing also a practical way to compute all these terms. This procedure is assessed on several test-cases computing the convergence curves obtained by using PC with respect to the reference solution, that is the exact analytical one.

Moreover, sensitivity indices based on skewness and kurtosis decomposition are introduced. The importance of ranking the predominant contributions in terms not only of the variance but also of higher order moments (thus extending the ANOVA analysis to high order statistic moments), is demonstrated with several test functions.

Two different strategies for reducing the complexity of the problem are considered: reducing the number of dimensions or limiting the order of interactions between different variables. Considering high-order statistics is shown to be of fundamental importance for saving the statistics properties of the reduced problem with respect to the complete one.

Future works will be directed towards computationally effective strategies for the reduction of the global computational cost, and the use of high-order statistics in robust design optimization following similar approaches already introduced for the variance [START_REF] Congedo | TSI metamodels-based multi-objective robust optimization[END_REF].

Proof. This term expresses the interaction between two multi-indices #» m p and #» m q , where one of them is raised to the second power. The two multi-indices should be different for construction, i.e. #» m p #» m q : if #» m q is not a subset of #» m p , then a set of coordinates #» m q \ ∩ pq belongs only to #» m q (if the set is totally disjointed, then the term ∩ pq is the null set). Note that the symbol ∩ pq indicates the coordinates contained in both the multi-indices p and q.

The integral reformulated into the form Ξpq

f 2 #» m p f #» m q dP pq = Ξp+∩pq f 2 #» m p        Ξq\∩pq f #» m q dP q\∩ pq        dP p+∩ pq = 0, (A.3)
where the internal integral is equal to zero due to the orthogonality of the ANOVA contributions (see equation ( 4)).

The case of the identification of the contributions to a specific multi-index #» m i is now addressed. The resulting is the set of the variables #» m pq = #» m i . If all the sub-sets of variables and their interactions of #» m i are collected in the set P i , the contributions to #» m i are a sub set of the 2 |P i | -1 simple combinations of class two. However, not all the combinations of sums #» m p ⊞ #» m q contribute to the conditional term µ 3 #» m i , for instance (1, 0, 0) ⊞ (1, 1, 0) #» m i . Thus, requiring simultaneously #» m q ⊂ #» m p and #» m pq = #» m i allows to identify the non-null contributions as

3 #» m p ∈P i #» m pq = #» m i P i ∋ #» m q ⊂ #» m p Ξpq f 2 #» m p f #» m q dP pq = 3 Ξi f 2 #» m i #» m q ∈P i, f #» m q dP i , (A.4)
where P i, is employed as shorthand for P i, = P i -{ #» m i }. The last term of (A.1) can be written as follows

6 N p=1 N q=p+1 N r=q+1 Ξpqr f #» m p f #» m q f #» m r dP pqr = 6 N p=1 N q=p+1 N r=q+1 #» m pq = #» m r Ξpq f #» m p f #» m q f #» m r dP pq . (A.5)
Proof. This case can be demonstrated extending what already done for the dyadic interaction between multi-indices:

Ξpqr f #» m p f #» m q f #» m r dP pqr = Ξpq+∩pqr f #» m p f #» m q       Ξr\∩pqr f #» m r dP r\∩ pqr       dP pq+∩ pqr = 0, (A.6)
by using the orthogonality property.

If a specific index #» m i is of interest, the conditional contribution is identified requiring #» m pqr = #» m i . Thus, the condition #» m pq = #» m r = #» m i identifies the non-null contributions

µ 3 #» m i = 6 #» m p #» m i #» m p #» m q #» m i #» m pq = #» m i Ξi f #» m p f #» m q f #» m i dP i . (A.7)
In conclusion, the final form for the third-order central moment µ 3 is equal to

µ 3 = N p=1 Ξp f 3 #» m p dP p + 3 #» m p #» m q ⊂ #» m p Ξpq f 2 #» m p f #» m q dP pq + 6 N p=1 N q=p+1 N r=q+1 #» m pq = #» m r Ξpq f #» m p f #» m q f #» m r dP pq , where each conditional contribution (µ 3 = µ 3 #» m i ) is µ 3 #» m i = Ξi f 3 #» m i dP i + 3 Ξi f 2 #» m i #» m q ∈P i, f #» m q dP i + 6 #» m p ∈P i, #» m p #» m q ∈P i, #» m pq = #» m i Ξi f #» m i f #» m p f #» m q dP i .
In the previous expression we note that the second and third term should include the product of #» m i and all the multiindices relative to the sub-interactions, i.e. P i . Moreover, P i, is used as shorthand for P i \ #» m i .

Appendix B. Fourth central moment expression for the ANOVA functional decomposition

In this section, Eq. ( 18) is obtained starting from the application of the multinomial theorem to Eq. ( 1)

µ 4 = N p=1 Ξp f 4 #» m p dP p + 4 N p=1 N q=1 q p Ξpq f 3 #» m p f #» m q dP pq + 6 N p=1 N q=p+1 Ξpq f 2 #» m p f 2 #» m q dP pq + 12 N p=1 N q=1 q p N r= j+1 r p Ξpqr f 2 #» m p f #» m q f #» m r dP pqr + 24 N p=1 N q=p+1 N r=q+1 N t=r+1 Ξpqrt f #» m p f #» m q f #» m r f #» m t dP pqrt . (B.1)
The first three terms are easy to handle. The first and the third one, on the right side, cannot be further simplified, while the second one can be manipulated as shown in the previous section. Thus, it is possible to write

4 N p=1 N q=1 q p Ξpq f 3 #» m p f #» m q dP pq = 4 #» m p #» m q ⊂ #» m p Ξpq f 3 #» m p f #» m q dP pq . (B.2)
as already demonstrated for the skewness term. More attention should be devoted to the last two terms. The first one of them can be written as

12 N p=1 N q=1 q p N r= j+1 r p Ξpqr f 2 #» m p f #» m q f #» m r dP pqr = 12 N p=1 N q=1 q p N r=q+1 #» m qr \∩ qr ⊆ #» m p Ξpqr f 2 #» m p f #» m q f #» m r dP pqr . (B.3)
Proof. If the multi-indices #» m q and #» m r are totally independent from the variables contained in #» m p , the condition Ξqr f #» m q f #» m r dP qr = 0 holds (due to the orthogonality of the ANOVA functional components). In the general case, when #» m qr ∩ #» m p 0, the existence of a null integral is related to the presence of variables in the multi-indices #» m q or #» m r which are not contained in #» m p . If #» m qr \ ∩ qr #» m p then it is possible to write as follows

Ξp +∩ qr f 2 #» m p        Ξ #» m qr\∩qr f #» m q f #» m r dP qr\∩ qr )        dP p = 0. (B.4)
Note that the internal integral is carried out with respect to a variable contained only in one of #» m q or #» m r , then is always zero due to orthogonality. Obviously, the case with the subsets related to #» m qr and #» m p totally disjointed, is included in the previous condition.

If a specific multi-index #» m i is provided, then in this case the contribution of this term is computed as follows 12

#» m p #» m p #» m q ∈P i #» m r ∈P i ,r>q #» m p ⊞∩ qr = #» m i Ξi f 2 #» m p f #» m q f #» m r dP i . (B.5)
Proof. The previous equation can be obtained considering the requirements #»

m pqr = #» m i and #» m qr \ ∩ qr ⊆ #» m p , with #» m p , #» m q , #» m r ∈ P i . It is easy to verify that if the second equation is true, then it must follow #» m p ⊞ #» m qr \ ∩ qr ⊆ #» m p ⊞ #» m p = #» m p , from which #» m i \ ∩ qr ⊆ #» m p . Finally, #» m i = #» m p ⊞ ∩ qr
holds (the equality sign follows from #» m p , #» m q and #» m r ∈ P i ). Remark that great attention must be paid in manipulating expressions with the summation of multi-indices ⊞. Generally, consider that #»

m p \ #» m q = #» m r ⇒ #» m p = #» m r ⊞ #» m q
holds but the contrary is not guaranteed.

The last term of the kurtosis can be also reformulated as follows

24 N p=1 N q=p+1 N r=q+1 N t=r+1 Ξpqrt f #» m p f #» m q f #» m r f #» m t dP pqrt = 24 N p=1 N q=p+1 N r=q+1 N t=r+1 #» m pq \∩ pq ⊆ #» m rt ⊆ #» m pq ⊞∩ rt Ξpqrt f #» m p f #» m q f #» m r f #» m t dP pqrt . (B.6)
Proof. This term can be obtained with some constraints on the multi-indices: they should share, two by two, some sets of coordinates: #» m pq \ ∩ pq ⊆ #» m rt and #» m rt \ ∩ rt ⊆ #» m pq . The fulfillment of the previous conditions assures that the integral cannot be divided into a product between integrals of orthogonal contributions. These conditions could be applied choosing randomly two couple of indices. Using the second constraint, #»

m rt ⊆ #» m pq ⊞ ∩ rt is obtained. Using this relation with the first requirement, #» m pq \ ∩ pq ⊆ #» m rt ⊆ #» m pq ⊞ ∩ rt holds.
If a set of variables is specified by using the multi-index #» m i , then the conditional contribution that arises from the previous term can be identified as follows 24

#» m p ∈P i #» m q ∈P i ,q>p #» m r ∈P i ,r>q t>r, #» m r ∈P i #» m i ⊆ #» m pq ⊞∩ rt #» m i ⊆ #» m rt ⊞∩ pq Ξi f #» m p f #» m q f #» m r f #» m t dP i . (B.7)
Proof. The set of all the possible sub-sets of variables relative to a multi-index #» m i is represented by P i , but not all the possible combinations of four elements selected from P i are leading to the multi-index #» m i . Thus, the first condition is to require that #» m pqrt = #» m i . In analogy with the previous proof, the non-null elements need to satisfy.

» m rt ⊞ #» m pq \ ∩ pq ⊆ #» m rt ⊞ #» m rt and #» m pq ⊞ #» m rt \ ∩ rt ⊆ #» m pq ⊞ #» m pq , the following conditions can be written #» m i \ ∩ pq ⊆ #» m rt ⇒ #» m i ⊆ #» m rt ⊞ ∩ pq #» m i \ ∩ rt ⊆ #» m pq ⇒ #» m i ⊆ #» m pq ⊞ ∩ rt . (B.8)
In conclusion, summing up all the contributions, the final form for the kurtosis can be written as where each conditional contribution, with respect to a fixed multi-index #» m i , is equal to

k = N p=1 Ξp f 4 #» m p dP p + 4 #» m p #» m q ⊂ #» m p Ξpq f 3 #» m p f #» m q dP pq + 6 N p=1 N q=p+1 Ξpq f 2 #» m p f 2 #» m q dP pq + 12 N p=1 N q=1 q p N r=q+1 #» m qr \∩ qr ⊆ #» m p Ξpqr f 2 #» m p f #» m q f #» m r dP pqr + 24
k #» m i = Ξi f 4 #» m i dP i + 4 Ξi f 3 #» m i #» m q ∈P i, f #» m q dP i + 6 #» m p ∈P i #» m p #» m q ∈P i #» m pq = #» m i Ξi f 2 #» m p f 2 #» m q dP i + 12 #» m p #» m p #» m q ∈P i #» m r ∈P i ,r>q #» m p ⊞∩ qr = #» m i Ξi f 2 #» m p f #» m q f #» m r dP i + 24 #» m p ∈P i #» m q ∈P i ,q>p #» m r ∈P i ,r>q t>r, #» m r ∈P i #» m i ⊆ #» m pq ⊞∩ rt #» m i ⊆ #» m rt ⊞∩ pq Ξi f #» m p f #» m q f #» m r f #» m t dP i .
(B.10)

Appendix C. Skewness from the PC expansion

In this section, the skewness is obtained from the PC series expansion. By applying the multinomial theorem, the skewness can be written as a sum of contributions generated by the interactions of the polynomial basis components At the first glance, it seems that no orthogonal contributions are present, because the interactions involve only polynomial forms raised to a power higher than one or triadic interaction. However, the second and third terms should be further investigated. Following from the definition of each polynomial term [START_REF] Novak | High dimensional integration of smooth functions over cubes[END_REF], the product between two polynomial terms of the basis, where the first one is raised to the power n with 1 < n ∈ N, can be written as follows

     dP = d i=1 Ξ i ψ n α p i (ξ i ) ψ α q i (ξ i ) p(ξ i )dξ i . (C.2)
Due to the orthogonality of the PC basis with respect to Ψ 0 = 1, it follows that if α p i = 0 then the integral with respect to the variable ξ i becomes

Ξ i ψ α q i (ξ i ) p(ξ i )dξ i = 0 for α q i 0. (C.3)
From this relation, the orthogonality of Ψ n p ( #» ξ ), Ψ q ( #» ξ ) follows. The non-null existence of the corresponding skewness (and kurtosis term) can be efficiently identified by means of the function ∆ p q defined in §3.2. The third term of the skewness from the PC series involves triadic interaction of polynomial terms raised to the power one

Ψ p ( #» ξ ), Ψ q ( #» ξ )Ψ r ( #» ξ ) = Ξ d         d i=1 ψ α p i (ξ i ) ψ α q i (ξ i ) ψ α r i (ξ i )         dP = d i=1 Ξ i ψ α p i (ξ i ) ψ α q i (ξ i ) ψ α r i (ξ i ) p(ξ i )dξ i . (C.4)
This term can be analyzed after the inspection of the relative multi-indices #» m ⋆,p , #» m ⋆,q and #» m ⋆,r . If the sum of the respective components of the multi-indices, i.e. m ⋆,p i + m ⋆,q i + m ⋆,r i , is equal to zero, then the variable is not present and no information can be obtained (the previous integral would be equal to 1 in such a case). If the sum m ⋆,p i + m ⋆,q i + m ⋆,r i is equal to 1, this means that the variable is present in only one polynomial term between ψ p , ψ q and ψ r , while it should not be present in the others (the relative coefficient α i = 0). This leads to a null integral due to the orthogonality of the basis with respect to the probability density function. However, another possibility can be associated to a null integral: if the sum m ⋆,p i + m ⋆,q i + m ⋆,r i = 2, the orthogonality between two polynomial terms guarantees that the integral is zero. The previous results can be resumed in the function ∆ pqr introduced in Section 3.2.

Appendix D. Kurtosis from the PC expansion

In this section, as already shown for the skewness in Appendix C, the kurtosis structure relying on the PC series is described. By applying the multinomial theorem to the PC series expansion, the kurtosis is computed as follows β p β q β r β t Ψ p Ψ q , Ψ r Ψ t .

(D.1)

The first and third term are not null in general, while the second one has already been analyzed for the skewness (the presence of the third power instead of the second one is not relevant). In addition, the last two terms need a separate analysis. The first contains the interaction between three polynomial terms, where the first of them is raised to the second power. In the general case of 1 < n ∈ N, it holds that

Ψ n p ( #» ξ ), Ψ q ( #» ξ )Ψ r ( #» ξ ) = Ξ d         d i=1 ψ n α p i (ξ i ) ψ α q i (ξ i ) ψ α r i (ξ i )         dP = d i=1 Ξ i ψ n α p i (ξ i ) ψ α q i (ξ i ) ψ α r i (ξ i ) p(ξ i )dξ i . (D.2)
From the last equation, it is possible to note that the orthogonality between polynomial term can be advocated if the term raised to the n-th power is zero. If α p i = 0 and the sum of the remaining terms m ⋆,q i + m ⋆,r i 0, then a null integral exist irrespectively of the coefficients α p i and α q i Ξ i ψ α q i (ξ i ) ψ α r i (ξ i ) p(ξ i )dξ i = 0. (D.3)

Note that the function ∆ p qr has been introduced in §3.3. The last term of the kurtosis expansion involves the interaction of four polynomials terms. This case represents an extension of the term already analyzed for the skewness (see Appendix C) where the interaction between three polynomial terms has been discussed. By inspecting the sum of the coefficients of the multi-indices #» m ⋆,p , #» m ⋆,q , #» m ⋆,r and #» m ⋆,t , it is possible to determine if the integral is equal to zero. In particular, if the sum relative to the ith coordinates, i.e. m ⋆,p i + m ⋆,q i + m ⋆,r i + m ⋆,t i , is equal to 1 or 2, the orthogonal properties with respect to the pdf holds (this is true irrespectively of the values of the α k i coefficients). This result has been used in Section 3.3 to define the function ∆ pqrt . The third and fourth-order central moments can be efficiently computed by means of the PC expansion if the functions ∆ q p , ∆ pqr , ∆ p qr and ∆ pqrt are evaluated before computing the integrals. Hereinafter, we refer to these functions as selecting functions. The first step is obtaining, in the pre-processing stage, the multi-indices #» α k ∈ N d for k = 0, . . . , P. This task can be accomplished using the algorithm proposed in [START_REF] Le Maître | Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics[END_REF]. Evaluating the selecting functions requires the computation of the normalized multi-indices #» m ⋆,k = #» m ⋆,k ( #» α k ). In the following we recalled the algorithm proposed in [START_REF] Le Maître | Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics[END_REF]

Appendix E. Computing the functions ∆

Figure 1 :

 1 Figure 1: Statistical moments error vs number of function evaluations in the case d = 2.

Figure 2 :Figure 3 :Figure 4 :Figure 5 :

 2345 Figure 2: Conditional statistics vs number of function evaluations in the case d = 2.

Figure 7 :

 7 Figure 7: Sensitivity indices for the first function f 1 (53) obtained with a PC series with total degree n 0 = 7.VariableTSI TSI s TSI k ξ 1 0.79 0.96 0.97 ξ 2 0.26 0.96 0.67 ξ 3 0.02 0.10 0.10

Figure 8 :

 8 Figure 8: SI for the function f 2 .

Figure 9 :

 9 Figure 9: PDFs for the complete g-function and the reduced models (see equations 54).

Figure 10 :

 10 Figure 10: PDFs for the complete g-function and the reduced models.

Figure 12 :

 12 Figure 12: PDFs for the complete f 2 and the reduced models up to the first and second orders.

Figure 13 :Figure 14 :

 1314 Figure 13: Computational grid for the turbine (a) and pressure contours for the reference case (b).

Figure 15 :

 15 Figure 15: Sensitivity indices for the Carnot factor (a) and the power output (b).

Figure 16 :

 16 Figure 16: Probability density functions for the metamodel including the first order interactions, second order interactions and the complete model. Efficiency is shown in (a) and the enthalpy variation in (b).

Figure 17 :

 17 Figure 17: Probability density functions for the metamodel including the first order interactions, second order interactions and the complete model. Carnot factor is shown in (a) and the power output in (b).

Figure 18 :Figure 19 :

 1819 Figure 18: Probability of exceeding a fixed value is computed as 1 -CDF for the enthalpy variation (a) and the power output (b).

Figure 20 :Figure 21 :Figure 22 :

 202122 Figure 20: Convergence of the constants C σ 2 , C s /(C σ 2 ) 3/2 and C k /(C σ 2 ) 2 varying the polynomial order n 0 for a metamodel including only first order interactions. The Carnot factor is reported in (a) while the power output is reported in (b).

  pq \∩ pq ⊆ #» m rt ⊆ #» m pq ⊞∩ rt Ξpqrt f #» m p f #» m q f #» m r f #» m t dP pqrt , (B.9)

  r Ψ 2 p , Ψ q Ψ r + 24

3 . Set P = d 4 . 5 .

 345 in which we add the computation of #» m ⋆,k :1. Set α 0 i = 0, for i = 1, . . . , d 1.bis #» m ⋆,0 = 0 2. Set α i j = δ i j , for 1 ≤ i, j ≤ d 2.bis #» m ⋆,i = #» α i Set h i (1) = 1 for i = 1, . . . , d For k = 2, . . . , n 0 a. Set L = P b. Set h i (k) = N ℓ=i h ℓ (k -1) c. For j = 1, . . . , d: * For ℓ = Lh j (k) + 1, . . . , L • Set P = P + 1 • Set α P i = α ℓ i for i = 1, . . . , d

Table 1 :

 1 Total sensitivity indices for the linear g-function function (52) based on a PC series with total degree n 0 = 5.

		TSI TSI s TSI k
	ξ 1	0.57 0.79 0.86
	ξ 2	0.29 0.56 0.64
	ξ 3	0.17 0.36 0.44
	ξ 4	0.11 0.24 0.31

Table 2 :

 2 Sensitivity indices for the quadratic g-function f gquad (52) obtained with a PC series with total degree n 0 = 5. Total sensitivity indices for the quadratic g-function f gquad (52) based on a PC series with total degree n 0 = 5.

		0.9
											variance
											skewness
		0.8									kurtosis
		0.7
		0.6
	Sensitivity index	0.4 0.5
		0.3
		0.2
		0.1
		0
		1	2	3	4	12	13	14	23	24	34 123 124 134 234 1234
								Conditional term
	Figure 6: Variable		TSI TSI s TSI k
						ξ 1		0.82 0.95 0.97
						ξ 2		0.14 0.47 0.44
						ξ 3		0.04 0.13 0.12
						ξ 4		0.01 0.04 0.04

Table 3 :

 3

Total sensitivity indices for the first function f 1 (53) based on a PC series with total degree n 0 = 7.

Table 4 :

 4 Total sensitivity indices for the first function f 2 (53) based on a PC series with total degree n 0 = 7.

		TSI TSI s TSI k
	ξ 1	0.36 0.70 0.71
	ξ 2	0.36 0.70 0.71
	ξ 3	0.36 0.70 0.71

Table 5 :

 5 Percentage

		Variance Skewness Kurtosis
	f G1	4.7997	29.236	15.039
	f G2	1.2369	7.7705	4.0632
	abs(µ n -µ n ex ) µ n ex

Table 6 :

 6 Total contribution for the variance, skewness and kurtosis up to the first and second order (total degree 5) as computed in Section 5.2 and their analytical counterparts.

	Function Variance Skewness Kurtosis
	f O1	86.46	8.02	7.83
	f O2	100.00	95.47	67.00
	f ex O1 f ex O2	82.76 98.78	0.00 81.32	31.51 76.52

Table 7 :

 7 Thermodynamic data for D 6 , where M is the percentage molecular weight, and T b is the boiling temperature at 1 atm. Properties are taken from Guardone et al.[START_REF] Guardone | Maximum intensity of rarefaction shock waves for dense gases[END_REF].

	M (g/mole) T c (K) P c (kPa) T b (K)
	444.9	645.8 961	518.1

Table 8 :

 8 Estimation errors for the linear metamodel following Section 4 for the turbine problem. The values relative to the amount of variance, skewness and kurtosis carried by the linear metamodel are reported in bold.

		35							0.0012
	Probability Density Function	5 10 15 20 25 30	Complete model First Order Second Order			Probability Density Function	0.0004 0.0006 0.0008 0.001 0.0002	Complete model First Order
										Second Order
		0
									0
		0.9	0.91	0.92	0.93	0.94	0.95	0.96	5000	5500	6000
					Efficiency				Dh
					(a)

In this case the ratios between the conditional central moment or the normalized one with the overall moment are equivalent.

We employ for it the Roman number I to avoid confusion with the first term of the ANOVA expansion.

Appendix A. Third central moment expression for the ANOVA functional decomposition

In this section, the skewness form Eq. ( 12) is computed using the expression obtained via the multinomial theorem applied to the functional ANOVA decomposition Eq. ( 1), and after the integration over Ξ it becomes

The previous equation displays interactions between variables and sub-sets of variables. In particular, the second and the third terms of could be simplified to highlight the contributions that are always equal to zero (due to the orthogonality).

The second term on the right hand side can be simplified as follows

Let us start with the evaluation of the function ∆ p q :

The selecting function ∆ pqr is evaluated as

The remaining functions, ∆ p qr and ∆ pqrt , are computed in analogy with the ones described above.