{'fraction_non_alphanumeric': 0.0303210092430677, 'fraction_numerical': 0.009118161378965775, 'mean_word_length': 5.110475100171723, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Common kernel-smoothed proper orthogonal decomposition (CKSPOD): An efficient reduced-order model for emulation of spatiotemporally evolving flow dynamics

In the present study, we propose a new surrogate model, called common kernel-smoothed proper orthogonal decomposition (CKSPOD), to efficiently emulate the spatiotemporal evolution of fluid flow dynamics. The proposed surrogate model integrates and extends recent developments in Gaussian process learning, high-fidelity simulations, projection-based model reduction, uncertainty quantification, and experimental design, rendering a systematic, multidisciplinary framework. The novelty of the CKSPOD emulation lies in the construction of a common Gram matrix, which results from the Hadamard product of Gram matrices of all observed design settings. The Gram matrix is a spatially averaged temporal correlation matrix and contains the temporal dynamics of the corresponding sampling point. The common Gram matrix synthesizes the temporal dynamics by transferring POD modes into spatial functions at each observed design setting, which remedies the phase-difference issue encountered in the kernel-smoothed POD (KSPOD) emulation, a recent fluid flow emulator proposed in Chang et al. (2020). The CKSPOD methodology is demonstrated through a model study of flow dynamics of swirl injectors with three design parameters. A total of 30 training design settings and 8 validation design settings are included. Both qualitative and quantitative results show that the CKSPOD emulation outperforms the KSPOD emulation for all validation cases, and is capable of capturing small-scale wave structures on the liquid-film surface faithfully. The turbulent kinetic energy prediction using CKSPOD reveals lower predictive uncertainty than KSPOD, thereby allowing for more accurate and precise flow predictions. The turnaround time of the CKSPOD emulation is about 5 orders of magnitude faster than the corresponding high-fidelity simulation, which enables an efficient and scalable framework for design exploration and optimization.

In the study of complex natural phenomena and engineering systems, high-fidelity simulations have been used for decades to provide detailed quantification of flow structures and dynamics that are otherwise hard to obtain from experiments or theoretical analyses [1][2][3][4][5][6][7]. These simulations , however, are computationally expensive and time consuming. For example, the simulation of propellant mixing in an azimuthal sector of a bi-swirl injector with 1.6 M cells takes 10,000 CPU hours (hexa-core AMD Opteron Processor 8431) to sweep 1 millisecond of flow time using largeeddy simulation (LES) [6]. It is thus impractical to rely solely on high-fidelity simulations for practical design, which often requires the survey of a wide parametric space. Surrogate models have been implemented to improve design efficiency. The objective of this study is to develop and validate an efficient surrogate model (emulator) to predict a spatiotemporally evolving flowfield in a faithful manner.
Surrogate models are designed to closely mimic the results of high-fidelity simulation with significantly reduced computing time and resource [8][9][10][11]. Surrogate models can be broadly categorized into data-fit, multifidelity (hierarchical), and reduced-order models (ROMs) [12].
Data-fit models are not physics-based and are formulated directly from interpolation or regression of simulation data. Such models fit a response surface connecting input and output data, by evaluating objective and constraint functions in the sampling space. Jones [13] presented a taxonomy of response surface models (RSMs) using a general function to represent different datafit models, including polynomial response surface models, kriging [14][15][16][17], and radial basis functions (RBFs) [18,19]. Another data-fit model is support vector regression [9], which can be viewed as an extension of RBFs; this model designates a threshold within which regression error in the sample data is acceptable, without affecting surrogate prediction. Data-fit surrogate models have been widely used in aerospace system design and optimization, but these models can face stiff challenges in problems with dynamic evolutions and high-dimensional data. Multifide lity (hierarchical) models can alleviate such an issue, by incorporating low-fidelity models into highfidelity models with a bridge function [20][21][22], where low-fidelity models are established via simplification of physical assumptions or reduction of numerical resolution. However, the prediction accuracy of multifidelity models may be compromised when the low-fidelity models provide a completely different trend from the high-fidelity counterpart.
ROMs are also frequently used to overcome the curse of dimensionality issue by constructing a low-dimensional subspace, on which the reduced operators or basis functions are sought from the original high-dimensional datasets. Projection-based ROMs are currently the most popular ROMs. This approach is often physics-based, since the reduced subspace typically contains the key structures of the dynamical system. Projection-based ROMs are particularly effective for systems whose input-output map is of low rank [23,24]. The most widely used projection-based ROM in engineering is the proper orthogonal decomposition (POD) [25]. POD is inherently connected to the idea of principal component analysis in the area of statistical learning [26] and Karhunen-Loève decomposition in the stochastic process literature [27]. In many problems, the POD-based ROM can capture key structures and dynamics embedded in high-fidelity data, which can then be used for prediction of a spatiotemporally evolving flowfield.
The use of POD-based ROMs has been discussed in many studies. In an early attempt, Ly et al. [28] applied POD to study the temperature field in a Rayleigh-Bénard convection problem, and used a cubic spline interpolation to predict the POD coefficients. Audouze et al. [29] employed radial basis functions to model POD coefficients, and validated the result against steady-state convection-diffusion-reaction problems. Swischuk et al. [30] proposed a physics-based parametric surrogate model using POD. In this paper, several machine learning methods, including neural networks, multivariate polynomial regression, k-nearest-neighbors, and decision trees, were used to learn the map between input parameters and POD expansion coefficients. They considered two engineering examples, and found that embedded physical constraints were important for the learned models. While these works have introduced some novel ideas for projection-based models, there are considerable limitations. First, these works considered only spatial or temporal development; no effort was given to problems with concurrent spatial and temporal evolution.
Further, all the methods proposed thus far have focused on training POD coefficients, with little attention to POD modes; this may be problematic when applied to cases with complicated dynamics and varying flow conditions and geometric parameters. The POD-based ROM has also been questioned for its intrinsically linear subspaces, although some nonlinearity may be embedded in POD expansion coefficients. Dynamic mode decomposition (DMD) [31,32] was developed to represent nonlinear, finite-dimensional dynamics without linearization by approximating the modes of the Koopman operator.
In recent years, deep learning methods, in particular, autoencoders, have been attempted to construct ROMs for engineering applications, due to their ability to treat system nonlinearity [33][34][35][36][37]. Autoencoders employ a neural-network structure with two elements: an encoder for the nonlinear mapping from the high-dimensional input to low-dimensional manifolds, and a decoder for nonlinear mapping from low-dimensional manifolds to an approximate representation of the high-dimensional input. Lee and Carlberg [35] implemented deep convolutional autoencoders to obtain reduced nonlinear manifolds of dynamical systems, which outperformed the linear subspace ROMs in selected advection-dominated problems. Xu and Duraisamy [37] proposed a three-level convolutional autoencoder network, including a convolutional autoencoder, a temporal convolutional autoencoder, and a fully-connected network, for parametric and future-state predictions of spatiotemporally evolving systems. Although deep learning methods embed nonlinear features through activation functions, the physical connection between reduced nonlinear manifolds trained by deep learning and actual dynamical structures remains unclear.
POD-based ROMs, on the other hand, provide a direct link between POD spatial modes and coherent structures of turbulent flows [25].
In this paper, we focus on the development of POD-based surrogate models, primarily due to the training data exhibits distinct physics. Every Gram matrix contains the temporal dynamics of the corresponding sampling point, and the common Gram matrix synthesizes these dynamics through element-wise multiplication. The CKSPOD technique is able to efficiently predict complex flowfields over a broad range of operating conditions and geometric parameters [40], and can be efficiently incorporated into a data-driven emulation framework for design analysis and optimization [41].
The present paper is structured as follows. Section II provides a detailed description of the CKSPOD methodology, including the proposed surrogate modeling framework and its training algorithm. Section III discusses the uncertainty quantification associated with the CKSPOD model.
In Section IV, the framework is applied for studying the spatiotemporal evolution of flow swirl injectors. Section V concludes this work with thoughts on future work.

The CKSPOD
In this section, we introduce a novel emulation technique, common kernel-smoothed proper orthogonal decomposition (CKSPOD). CKSPOD provides an efficient way to train a projectionbased surrogate model using simulation results at observed design settings, and allows for predictions (emulations) over the desired design space in practical turn-around times. It also circumvents the issue of phase differences (sign differences of eigenvectors) commonly found in the KSPOD method by constructing a transfer matrix in the data reduction process.
The fundamentals of POD are briefly introduced first. POD is a model reduction method that extracts orthogonal basis functions and associated temporal coefficients based on the energy norms. (2) Here, is an × orthonormal matrix spanning the column space of , an × diagonal matrix of singular values, and an × orthonormal matrix spanning the row space of . and can be computed from the eigenvalues and eigenvectors of and , respectively. Let be the inner product of the data matrix , also known as the Gram matrix.
The eigenvectors of make up the columns of , i.e., = , where = ∈ × is the diagonal matrix of eigenvalues of .
All POD modes and time-varying coefficients of can then be written as: = , Model reduction is possible when the number of POD modes selected for reconstruction is truncated at a lower rank of ( < ). The POD technique is frequently used to develop projection-based surrogate models in computational fluid dynamics.
As mentioned in the Introduction, CPOD-and KSPOD surrogate models have been proposed for problems with varying flow conditions and geometries [17,38]. CPOD concatenates all training information into a large matrix, = ( , = 1, 2, ⋯ , ) ∈ × , which is then used to build the CPOD covariance matrix and CPOD modes. To do this, a physics-guided common grid system is designed for the projection of the original database. A set of CPOD modes are obtained following SVD, and the corresponding time-varying coefficients at a new design setting are determined by the kriging procedure [17,38]. On the other hand, the dimension of these CPOD modes ( ) is significantly larger than that of POD modes ( ) at each training setting, producing an averaging effect on the flowfield. Consequently, the prediction accuracy of instantaneous flow information is compromised.
To overcome the limitations of CPOD, the KSPOD technique constructs spatial functions using the kriging-weighted average of POD modes at every sampling point, retaining dominant structures across all design points [39]. The method considerably improves prediction of instantaneous flowfield. A major difficulty, however, arises when the POD modes at different sampling points deviate in phase, as manifested by the element signs in the POD mode matrix. To this end, we propose a CKSPOD emulation to circumvent the issue of phase shift and provide faithful prediction of spatiotemporal flow dynamics in the present study. We now discuss the details of the proposed methodology.
In CKSPOD, a common Gram matrix ℂ, which contains data information of all training cases, is constructed as the Hadamard product of the Gram matrices of these cases. The Hadamard operator, denoted as ∘ , is the element-to-element multiplication of two matrices of similar dimension, as follows.
Accordingly, the common Gram matrix can be written as where is the column matrix of eigenvectors of ℂ and the corresponding diagonal matrix. The Gram matrix , the spatially averaged temporal correlation matrix, contains the temporal dynamics of the corresponding sampling point , and the common Gram matrix synthesizes these dynamics through element-wise multiplication using the Hadamard product. This operation can adjust the phase-shift issue encountered in the previously proposed KSPOD model. Note that an inherent assumption applied in Eq. (7) is that the number of snapshots (m) collected for each design setting is identical. This ensures that all Gram matrices are of the same dimension, and enables the Hadamard product. For problems with different geometries and computational grids, an additional procedure is needed to project the original simulation data into a common grid space [17].
If we also define Π as the notation of the Hadamard product, ℂ can be organized as: Directly using the POD modes from the original Gram matrix of each training case, as in KSPOD, the CKSPOD technique implements the common Gram matrix in Eq. (8) to deduce spatial functions and time-varying coefficients for the training cases. The transferred spatial functions of the design setting i, denoted as ′ , are written as Substituting ℂ using Eq. (8), ′ can be related to the original POD modes in the following Here the matrix is defined as We call the CKSPOD transfer matrix, which converts the original POD modes to the spatial functions at each design setting. To utilize these spatial functions properly with kriging-based weights for the new design setting, a normalization process is performed.
Accordingly, the transferred matrix of time-varying coefficients is expressed as Once the transferred matrices of new spatial functions and time-varying coefficients are deduced, kriging is implemented to develop the basis functions and coefficients for a new design setting. It should be noted that the transferred spatial functions ̃ maintain the orthogonality after the conversion by the CKSPOD transfer matrix, and can thus be considered as new POD modes.

The kriging model
Kriging, also known as Gaussian process (GP) regression, models the responses over the input space as a sample drawn from a GP. For a set of design settings { ∈ } =1 , the observed functions of interest are weighting parameters of spatial basis modes and corresponding mode coefficients. Based on the training dataset of the − pair, with a new input , kriging predicts the corresponding response . The mathematical formula of kriging in terms of prediction is given by where ̂= − / − is the maximum likelihood estimate of , is a -vector of 1's, and is a × matrix of re-parameterized squared-exponential correlation function whose vector whose i-th entry is ( , ). This allows for a more numerically stable optimization of maximum likelihood estimators [17]. The spatial basis functions are calculated in a slightly different way. Kriging is used to predict the weight of each spatial function at observed points on the spatial function at the new design setting.
To this end, the observations are now taken to be the orthonormal vector , where is avector with 1 in its -th element and 0 elsewhere. Intuitively, this quantifies the fact that the spatial mode information extracted in the -th design setting corresponds to only that setting and not the other − 1 settings. With this in mind, the resulting predictor in Eq. (14) can be viewed as the predicted weight for that particular spatial mode at the new design setting , denoted as ̂, This procedure is repeated for each of the unit vectors ( ) =1 , from which the weighting parameters (̂, ) =1 can be obtained. The weighting parameters are normalized to ensure their summation equal to unity.
The weighting parameters are subsequently used to predict the new spatial function modes through a weighted average of the extracted modes at the new design settings, expressed as, With the time-varying coefficients and spatial functions obtained using Eqs. Note that the KSPOD approach [39] directly uses the original POD modes ( ) and coefficient ( ) at each design setting i to perform the kriging process using Eq. (14).

TRAINING:
Step 1: Construct the common Gram matrix ℂ based on Gram matrices of all observed design settings using Hadamard product (Eq. (8)).
Step 4: Perform the kriging model with inputs { 1 , … , } and the predictive function for is the weight ̂( ) associated with each design setting i. The predicted spatial function of mode k, ̂( , ) is the weighted sum of corresponding spatial functions at observed design settings (Eq. (15)).

PREDICTION:
Use Eq. (16) to construct the spatiotemporal flowfield data at the new design setting.

Uncertainty quantification
For surrogate modeling, it is also crucial to quantify the uncertainty of the prediction to assess the model accuracy. Here, we assume that the database created by high-fidelity LES simulations is reliable and accurate, and the uncertainty of statistical prediction primarily results from the kriging process. Our previous work in CPOD has shown that invoking the conditional distribution of the multivariate normal distribution, the kriging-predicted time-varying coefficients at a new design setting follow the Gaussian distribution, Here the minimum mean square error (MMSE) predictor ̂ for ( )|{̃( )} =1 and its corresponding variance are given by, where , , and are × identity matrix, 1-vector of q elements, and the × covariance matrix, respectively.
Similarly, the variance associated with the weights of the predicted spatial functions during kriging at a new design setting can also be represented by Eq. (21), except with as an identity matrix.
The UQ of the final prediction using Eq. (16)  been extensively studied [43,44], and these methods are utilized for computing the pointwise CI of ( , ) in Section III. With this in hand, the prediction and UQ of TKE from the CKSPOD and KSPOD are compared with the simulated TKE at the validation design point.

Model validation: Swirl injectors
The proposed CKSPOD emulation methodology is applied to the design of a swirl injector.
Swirl injectors have been widely used to achieve efficient mixing and combustion in propulsion and power-generation systems [1,16]. The performance of the CKSPOD technique is examined in detail in this section. Figure 1 shows sectional views of the simplex swirl injector of concern [16][17][18], identical to the one used in our previous KSPOD work [39]. Liquid oxygen (LOX) is tangentially introduced into the injector through orifices. In this study, the orifices are replaced by a slit on the injector wall;

Swirl injector configuration
the slit width ( ) is carefully chosen to ensure a mass flow rate identical to that of the discrete orifices. A liquid film forms along the injector wall [45][46][47][48]. A hollow gaseous core is formed in the center region due to the conservation of angular momentum. The liquid film exits the injector as a thin conical sheet and subsequently undergoes rapid atomization. The detailed flow dynamics of this type of simplex swirl injector under supercritical conditions has previously been thoroughly studied using LES [46,47].

Fig. 1. Schematic of swirl injector
The selection of injector design variables requires careful assessment of their influences on flow dynamics. The major design attributes include the injector length, L, radius, R n , tangential inlet width, δ, tangential injection angle, θ, and distance between the inlet and headend, ΔL. A sensitivity analysis described in previous studies [17,38] has identified δ, θ, and ΔL as the most significant parameters. They are thus selected as the design parameters in the current study. The baseline geometry and operating conditions of the injector, including the LOX inlet temperature T in , ambient temperature T ∞ , ambient pressure p ∞ , and mass flow rate ṁ, are tabulated in Table 1.
The range for each design variable is listed in Table 2. The distance between the inlet and headend, ΔL, is typically 1.5-2.0 times the diameter of the injection orifice [45]; the location depends on a tradeoff to avoid (1) excessive viscous losses when the injection slit is too close to the headend, and (2)  With the selected design parameters and their ranges, the SLHD method is implemented to generate the design settings with the =10p rule-of-thumb for sample size [49]. A total number of 30 design settings are employed in the present study, as listed in In SLHD, the space-filling property of the design points in each slice is optimal. The overall design matrix contains five slices, and each slice includes six design settings. Fig. 2 shows twodimensional projections of the design settings categorized by different slices. LES-based highfidelity simulations are performed at all design settings.

Fig. 2. Two-dimensional projection of design points obtained by Sliced Latin Hypercube
Design (SLHD) methodology in the design space.

High-fidelity simulation
The theoretical formulation for high-fidelity simulations is described in detail in [50]. It allows for detailed modeling of mixing and combustion over the entire range of fluid thermodynamic states of concern [51]. Turbulence closure is achieved by means of an LES technique.
Thermodynamic properties are evaluated according to fundamental thermodynamics theories and a modified Soave-Redlich-Kwong (SRK) equation of state (EOS). Transport properties, including thermal conductivity and viscosity, are determined using extended corresponding-state principles.
Mass diffusivity is estimated based on the Takahashi method, calibrated for high-pressure conditions.
The numerical framework is based on a preconditioning scheme with a unified treatment of general-fluid thermodynamics [50,52]. It utilizes a density-based, finite-volume methodology, along with a dual-time-step integration technique. Temporal discretization is accomplished by a second-order backward difference; and the inner-loop pseudo-time term is integrated via a fourstep Runge-Kutta scheme. A fourth-order central difference scheme in generalized coordinates is used to obtain spatial discretization. Fourth-order matrix dissipation is applied to secure numerical stability with minimum contamination of the solution. Finally, a multi-block domain decomposition technique associated with the message passing interface technique for parallel computing is applied to optimize computation performance.

POD modes and CKSPOD spatial functions
With spatiotemporal simulation data assigned for all 30 design settings, a database is established to train the emulator. To justify the applicability of CKSPOD, the original POD modes are matched to the transferred spatial functions deduced from the common covariance matrix, aliased as the CKSPOD modes for convenience. Figure 3 shows the first four modes of POD and CKSPOD for the pressure field of Case 16. All modes are normalized to facilitate comparison. The POD and CKSPOD modes exhibit similar features, which suggests that most significant coherent flow structures in the original design setting are retained after the CKSPOD procedure. The energy percentage of each POD mode, however, is consistently larger than its CKSPOD counterpart.
Mathematically, the energy percentage is represented by the ratio of the diagonal element (eigenvalue) to the trace of the diagonal matrix . The energy carried by the fist POD mode is around 2.5 times of the corresponding CKSPOD mode; this can be explained by the introduction of the CKSPOD transfer matrix defined in Eq. (11). This transfer matrix contains the dominant dynamics embedded in the other 29 design settings, and therefore reduces the energy percentage of the original POD mode.

modes.
To further demonstrate the similarity, Figure   Fig. 4 shows probability density distributions of the first four pressure modes of POD and CKSPOD in Case 16 based on the kernel smoothing function. A kernel distribution is a nonparametric representation of the probability density function, ℎ ( ), of a random variable, written as where the sample size, the kernel density smoothing function, and the bandwidth acting as a smoothing parameter. Here, the bandwidth of the kernel smoothing function is optimally selected for estimating densities for the normal distribution [53] to produce a reasonably smooth curve. Close similarity is observed between the probability density distributions of the POD and CKSPOD modes. Both Modes 1 and 3 show negative skew with a tail on the left. The distributions of the fourth POD and CKSPOD modes are almost overlapped, with the latter slightly shifted to the left. The transfer matrix changes the mean value of the distribution and rescales the magnitude the modes. Such behavior is critical so POD modes for all design settings can be automatically adjusted, unlike the situation with the KSPOD approach [39], in which all POD modes must be manually checked to avoid the occurrence of the kriging procedure.

Fig. 4. Probability densities of POD Modes 1-4 for Case 16 from Cluster A and corresponding modes. Vertical lines represent mean values and represents bandwidth of kernel
smoothing function. Figure 5 shows the accumulated energy percentage of the POD modes for the pressure field for 8 randomly selected cases from the 30 design settings. Also included are the results for the CKSPOD modes. In the CKSPOD approach, all of the 30 cases are treated using the identical eigen-decomposition process with the same energy accumulation through the CKSPOD transfer matrix. The phase-difference issue encountered in the KSPOD method is thus remedied. On the other hand, the energy accumulation is much slower for the CKSPOD modes, due to the data smoothing effect. It takes the first 173 modes to obtain 90% of the total energy for the CKSPOD method, compared to the first 23-99 modes for the POD method. More modes are required for the CKSPOD approach to capture the same amount of energy in the establishment of the emulator. In the present study, the first 522 modes, which cover more than 99% percentage of the total energy, are included to the train the surrogate model and build the spatial functions and coefficients at the new design setting.  Table 4. Each group consists of two validation cases to fully evaluate the performance of the CKSPOD emulation. High-fidelity simulations are also performed at these settings to validate the predicted results by the emulation.  Figure 6 shows the instantaneous density fields of the LES-based simulation and the CKSPODemulation for Case D2. The evolution of the liquid film (seen in the dense fluid) and its spreading downstream of the injection port agree well between the simulation and emulation.

Temporal evolution
The performance of the CKSPOD method is further examined in terms of the temporal evolution of the flowfield. Figures 11 and 12 show the density field at three different time instants for the KSPOD and CKSPOD emulations, respectively, for Case A2. Also included are the LES predictions. The KSPOD method is able to predict the rolling vortices upstream of the injector ( ≤ 15 mm), but produces weakened surface structures. The phenomenon is attributed to the possible difference among the POD modes of the training dataset, which leads to the cancellation of flow dynamics in kriging and reconstruction [12]. The CKSPOD emulation, shown in Fig. 12, on the other hand, can capture most significant rolling vortices and the stringy ligaments revealed in the LES results at all time instants. The traveling surface wave propagates downstream all the way to the injector exit at the speed predicted by the LES method.
To establish broad confidence, the validity of the CKSPOD method is assessed against all test cases. Fig. 13 shows

Performance measures
Two performance metrics, liquid film thickness and spreading angle at the injector exit, are calculated to quantitatively assess the model accuracy. Table 5 lists the time-mean values, standard deviations, and relative errors obtained from the LES, KSPOD, and CKSPOD methods for all validation cases. Time-mean values are determined by averaging the instantaneous data over a statistically meaningful duration. SD denotes standard deviation, and the error is calculated as follows.
where represents data from simulation and from emulation. Both the KSPOD and CKSPOD emulators are capable of producing high-fidelity results of time-mean flow quantities.
The CKSPOD method, however, outperforms the KSPOD method, with errors on the order of 0.1% for most cases. The spatial distributions of time-mean liquid film surface along the axial direction are also obtained by averaging more than 1000 snapshots. Figure 14 shows the results along with the LES data. The CKSPOD-predicted liquid film surface closely coincides with the LES predictions, while the KSPOD results show discrepancies near the transient region in Cases B1, C1, and D1.

Uncertainty quantification (UQ)
For computer experiments, quantification of prediction uncertainty is as important as the prediction itself. In our earlier study, the one-sided width of the 80% pointwise confidence interval of axial velocity and turbulent kinetic energy was used as a measure for the spatial uncertainty quantification [38]. The work demonstrated the usefulness of uncertainty quantification not only as a measure of predictive uncertainty but also as a means for extracting useful flow physics without expensive simulations. In this subsection, the uncertainty quantification of both KSPOD and CKSPOD emulations are explored. The turbulence kinetic energy defined in Sec. 3 is analyzed.  [1,4]. The turbulent kinetic energy predicted by the CKSPOD method bears close resemblance to that predicted by the LES method. A very similar recirculating flow in the center is shown between CKSPOD and LES cases. The KSPOD emulation leads to considerable over-prediction for Cases A1 and C1. The situation is further corroborated by the standard deviation, as shown in Fig. 17.
The smaller the standard deviation, the lower the uncertainties. The CKSPOD method outperforms KSPOD, with a much smaller standard deviation in turbulent kinetic energy for both validation cases. Similar observations are made for all other cases and flow quantities (not shown here).  Therefore, the time saving of emulation is more than 5 orders of magnitude, compared to simulation. The difference between the KSPOD and CKSPOD methods lies in the execution of eigendecomposition; the later takes a slightly more CPU time due to extra procedures for building the common Gram matrix and transfer functions. The overall wall time for the CKSPOD emulator is about 50 seconds of CPU time per snapshot, roughly 1.2 times longer than the KSPOD emulation.