## 1 Introduction  
Over the past few decades, modern computational methods have been very successful in predicting the evolution of systems of great complexity. However, in real-life applications, the system behavior is significantly affected by parametric uncertainties. As an example, geometric uncertainties in the shape of the front wing, diffuser, and the tip of a Formula (1) can affect its performance [1]. A small selection of other examples can be found in References [2; 3; 4; 5]. Various techniques for the quantification of uncertainties have been developed. These include sampling techniques, such as the Monte-Carlo [6], Latin Hypercube sampling [7], and quasi Monte-Carlo [8]. However, despite their accuracy and ease of application, these methods become computationally intractable for large, real world problems, such us high fidelity simulations of turbulent flows. Other approaches, based on the spectral representation of uncertain quantities, have also been developed. These techniques are referred to as Polynomial Chaos Expansion (PCE) methods. They were first developed in Reference [9] for Gaussian uncertainties, and, later in Reference [10], they were generalized for all Probability Density Functions (PDFs) in the Askey scheme. The ability of generalized PCE (gPC) in quantifying uncertainties has been amply demonstrated in many applications, for example, in Computational Fluid Dynamic (CFD) systems [11; 12; 13]. Recent developments in Uncertainty Quantification (UQ) include using the PCE in data-driven approaches [14]. Spectral approaches are computationally efficient for a small number of uncertain parameters, \(m\), but their cost increases exponentially with \(m\) (curse of dimensionality). There are two major approaches to applying gPC, the intrusive (iPC) [5; 15] and the non-intrusive (niPC) [16]. In the intrusive approach, the state vector is written in spectral form using the gPC, and a new set of state equations, the iPC equations, are extracted using Galerkin projection. On the other hand, the non-intrusive approach is based on the spectral representation of a Quantity of Interest (QoI), such as the lift of an airfoil, and the spectral coefficients are calculated by repeated evaluations of the QoI at specific quadrature locations. For general unsteady systems, the iPC approach is not feasible for conducting UQ because the accuracy of the spectral representation of the state variables starts to deteriorate after the system starts to evolve from the initial condition. For this reason, we used the niPC approach and focused on uncertainties of time-averaged quantities and their sensitivities with respect to input parameters. Time-averaged quantities are of interest to many engineering applications, for example, in aeronautics (lift and drag around an airfoil), energy generation (fuel consumption and pollutant formation), process industries (average heat transfer, mixing rate), etc. Although time-averages of state variables can be easily computed by collecting statistics over long time integration, the evaluation of their sensitivities to input parameters is much more involved. The adjoint sensitivity analysis approach [17] is based on a linearization of the evolution equation with respect to one or more system parameters. This approach, however, breaks down in chaotic systems, where the state variables display extreme sensitivity to minute variations in input conditions [18]. In such cases, the adjoint variables diverge; this is popularly known as the 'butterfly effect', and it is due to the presence of one or positive Lyapunov exponents. The Least-Squares Shadowing (LSS) approach was developed in Reference [19] and resolves this problem. It is based on the shadowing lemma for ergodic and uniformly hyperbolic dynamical systems [20]. Although computationally expensive, the method can give values for the sensitivities of time-averaged quantities that match well with finite differences.