## 2 Materials and Methods  
### Ordinary Kriging  
Kriging was first introduced by [PERSON] et al. as a surrogate for expensive computational models [26]. Given \(n\) initial training points, the response of the OK model can typically be expressed as follows:  
\[\hat{y}(x)=\hat{\mu}+\hat{e}(x) \tag{1}\]  
Here, \(\hat{\mu}\) represents the mean term, which is determined by solving a weighted least squares problem:  
\[\hat{\mu}=\left(F^{\mathrm{T}}R^{-1}\mathrm{F}\right)^{-1}\mathrm{F}^{\mathrm{ T}}R^{-1}\mathbf{y} \tag{2}\]  
Here, \(\mathbf{F}\) is the regression matrix representing the zero-order basis function. In the OK model, the polynomial vector of the basis function can be expressed as follows:  
\[\mathbf{F}=\mathbf{1}=\underbrace{\left[\mathbf{1},\cdots,\mathbf{1}\right]}_{n}^{\mathrm{T}} \tag{3}\]  
The correlation matrix \(R\) characterizes the spatial correlation among all training points:  
\[R=\begin{bmatrix}\psi\left(\mathbf{x}^{1},\mathbf{x}^{1}\right)&\ldots&\psi\left(\mathbf{ x}^{1},\mathbf{x}^{n}\right)\\ \vdots&\ddots&\vdots\\ \psi\left(\mathbf{x}^{n},\mathbf{x}^{1}\right)&\ldots&\psi\left(\mathbf{x}^{n},\mathbf{x}^{n} \right)\end{bmatrix} \tag{4}\]  
Here, the Gaussian function is selected as the modeling function for spatial correlation:  
\[\psi\left(x^{i},x^{j}\right)=exp\left[\sum_{k=1}^{n}-\theta_{k}\left(\left\| x_{k}^{i}-x_{k}^{j}\right\|\right)^{2}\right] \tag{5}\]  
\(\theta_{k}\) is a critical hyper-parameter in kriging, and its value directly influences the quality of the kriging fit. As seen from the expression, \(\theta_{k}\) scales the influence of the \(k_{th}\) design variable on the correlation between training points. Figure 2 illustrates the variation in correlation with respect to the distance between training points under different hyper-parameter values. The determination of hyper-parameters requires optimization using the maximum likelihood estimation (MLE) method to achieve the optimal fitting trend for the training points.  
The error term \(\hat{e}(x)\) in Equation (1) is treated as a realization of a stochastic process and can be expressed as follows:  
\[\hat{e}(x)=\mathbf{w}^{T}\mathbf{r}(x) \tag{6}\]  
Here, the correlation vector \(\mathbf{r}\) describes the spatial correlation between the training points and the prediction point:  
\[\mathbf{r}=\left[\psi\left(\mathbf{x}^{n+1},\mathbf{x}^{1}\right),\ldots,\psi\left(\mathbf{x} ^{n+1},\mathbf{x}^{n}\right)\right] \tag{7}\]  
The determination of \(\mathbf{w}\) requires solving a system of linear equations:  
\[\mathbf{w}=\mathbf{R}^{-1}(\mathbf{y}-\mathbf{F}\hat{\mathbf{\mu}}) \tag{8}\]  
The final kriging predictor can be expressed as follows:  
\[\hat{\mathbf{y}}\left(\mathbf{x}^{n+1}\right)=\hat{\mathbf{\mu}}+\mathbf{w}^{T}\mathbf{r} \tag{9}\]  
The mean square error of the prediction made by the kriging predictor is given by the following:  
\[s^{2}\left(x^{n+1}\right)=\hat{\sigma}^{2}\left[1-\mathbf{r}^{T}\mathbf{R}^{-1}\mathbf{r} +\frac{\left(1-\mathbf{r}^{T}\mathbf{R}^{-1}\mathbf{r}\right)^{2}}{\mathbf{r}^{T}\mathbf{R}^{-1} \mathbf{F}}\right] \tag{10}\]  
Here, \(\hat{\sigma}^{2}\) represents the predicted system variance:  
\[\hat{\sigma}^{2}=\frac{1}{n}(\mathbf{y}-\mathbf{F}\mu)^{T}\mathbf{R}^{-1}(\mathbf{y}-\mathbf{F}\mu) \tag{11}\]  
The objective function values at the unknown points can be predicted by solving the linear equation systems (2), (8), (10), and (11).