## Appendix A Polynomial Chaos Expansion
#### Treating Time Dependency  
Suppose that the computational model \(\mathcal{M}\) reports the response \(\mathbf{y}=\{y_{1},...,y_{n_{out}}\}\) for a given \(\mathbf{\theta}\) at a fixed set of time instances \(\mathbf{T}=\{t_{1},...,t_{m}\}\). Since the cost of training and prediction using a surrogate model is \(\mathcal{O}\big{(}(n_{t}n_{out})^{3}\big{)}\), even moderate values of \(n_{out}\) and \(n_{t}\) renders the computations difficult. Also, individual surrogate modeling may result in a high degree of redundancy, since the simulation outputs in the form of time series are highly correlated [66]. To circumvent this computational issue, one can treat the temporal dependency of the outputs via a dimensionality reduction technique that was first introduced by Reference [67]. We employ Principal Component Analysis (PCA) to obtain a linearly uncorrelated variables with decreasing variances by applying an orthogonal transformation.  
Let us consider a random vector \(\mathbf{\tilde{Y}}\in\mathbb{R}^{n_{t}\times n_{out}}\) containing the \(Q\) realization of model outputs \(\mathbf{Y}=\{\mathbf{y}^{(1)},...,\mathbf{y}^{(Q)}\}^{\top}\). The exact mean ( \(\mathbf{\mu}_{\mathbf{\tilde{Y}}}=\mathbb{E}[\mathbf{\tilde{Y}}]\)) and the covariance (\(\mathbf{\Sigma}_{\mathbf{\tilde{Y}}}=\text{Cov}[\mathbf{\tilde{Y}}]\)) are approximated by:  
\[\bar{\mathbf{\mu}}_{\mathbf{\tilde{Y}}}=\frac{1}{Q}\sum_{q=1}^{Q}\mathbf{y}^{(q)}, \quad\mathbf{\Sigma}_{\mathbf{\tilde{Y}}}=\frac{1}{Q-1}\sum_{q=1}^{Q}(\mathbf{y}^{(q) }-\bar{\mathbf{\mu}}_{\mathbf{\tilde{Y}}})(\mathbf{y}^{(q)}-\bar{\mathbf{\mu}}_{\mathbf{ \tilde{Y}}})^{\top}.\] (A14)  
For \(q=1,...,Q\), the PCA transformation can be performed by:  
\[\mathbf{\tilde{\xi}}^{(q)}=\mathbf{\Phi}_{N^{*}}^{\top}(\mathbf{y}^{(q)}-\bar{\mathbf{\mu}}_{ \mathbf{\tilde{Y}}}),\] (A15)  
where \(\mathbf{\bar{\Phi}}_{N^{*}}^{\top}=(\bar{\Phi}_{0},...,\bar{\Phi}_{N^{*}})\) is composed of \(N^{\star}\) eigenvectors for which the proportion \(\Sigma_{i=1}^{N^{*}}\bar{\lambda}_{i}/\Sigma_{i=1}^{N^{*}}\bar{\lambda}_{i}\) of the total empirical variance is equal or larger than a pre-specified threshold. The eigenvectors \(\bar{\Phi}_{i}\) and the eigenvalues \(\bar{\lambda}_{i}\) must fulfill \(\mathbf{\Sigma}_{\mathbf{\tilde{Y}}}\bar{\Phi}_{i}=\bar{\lambda}_{i}\bar{\Phi}_{i}\) for \(i=1,...,N\). The reduced model output can be represented by \(\mathcal{Z}=\{\mathbf{\tilde{\xi}}^{(1)},...,\mathbf{\tilde{\xi}}^{(Q)}\}^{\top}\) and reconstruction of the originally observed samples for \(q=1,...,Q\) can be approximated as:  
\[\mathbf{y}^{(q)}\approx\bar{\mathbf{\mu}}_{\mathbf{\tilde{Y}}}+\mathbf{\bar{\Phi}}_{N^{*} }^{\top}\mathbf{\tilde{\xi}}^{(q)}=\bar{\mathbf{\mu}}_{\mathbf{\tilde{Y}}}+\sum_{i=0} ^{N^{*}}\xi_{i}^{(q)}\Phi_{i},\] (A16)  
The compressed data set \(\mathcal{Z}=\{\mathbf{\tilde{\xi}}^{(1)},...,\mathbf{\tilde{\xi}}^{(Q)}\}\) for only 200 realization of the thermochemical model, introduced earlier, are visualized in the parallel coordinate plot in Figure A1, only for illustration purposes. It is evident that the main components are ordered according to their individual contribution to the total variance and centered at zero.  
Once the model outputs are transformed, we only need to train the surrogate model to map the unknown parameters \(\mathbf{\theta}\) with the components of the forward model responses, as shown in Figure A1. Then, Equation (A16) is used to obtain the prediction of the original model responses for the task of inference.