## 2 p-Refined Multilevel Quasi-Monte Carlo
### Modeling the Spatial Variability  
The uncertain spatial variability of the Young's modulus (model problem 1) and the soil's cohesion (model problem 2) is represented by means of a random field. The Young's modulus is represented as a gamma random field, while the soil's cohesion is represented as a lognormal random field. The starting point for the construction of both fields is the generation of a (truncated) Gaussian random field. This is done by using a Karhunen-Loeve (KL) expansion, see [32]. Consider a Gaussian random field \(Z(\mathbf{x},\omega)\), where \(\omega\) is a random variable with a given covariance kernel. For both cases we choose the Matern covariance kernel,  
\[C(\mathbf{x},\mathbf{y}):=\sigma^{2}\frac{1}{2^{\
u-1}\Gamma\left(\
u\right)} \left(\sqrt{2\
u}\frac{\|\mathbf{x}-\mathbf{y}\|_{2}}{\lambda}\right)^{\
u}K _{\
u}\left(\sqrt{2\
u}\frac{\|\mathbf{x}-\mathbf{y}\|_{2}}{\lambda}\right)\,, \tag{21}\]  
with \(\
u\) the smoothness parameter, \(K_{\
u}\) the modified Bessel function of the second kind, \(\sigma^{2}\) the variance and \(\lambda\) the correlation length. The parameter values, \(\
u\), \(\lambda\) and \(\sigma\) will be given later on for both cases. The KL expansion can be formulated as:  
\[Z(\mathbf{x},\omega)=\overline{Z}(\mathbf{x})+\sum_{n=1}^{\infty}\sqrt{\theta _{n}}\xi_{n}(\omega)b_{n}(\mathbf{x})\,. \tag{22}\]  
Here, \(\overline{Z}(\mathbf{x})\) denotes the mean of the field, and is set to zero. The \(\xi_{n}(\omega)\) denote i.i.d. standard normal random variables. The symbols \(\theta_{n}\) and \(b_{n}(\mathbf{x})\) respectively denote the eigenvalues and eigenfunctions of the covariance kernel corresponding to Equation (21). They are the solutions of the eigenvalue problem  
\[\int_{D}C(\mathbf{x},\mathbf{y})b_{n}(\mathbf{y})\mathrm{d}\mathbf{y}=\theta_ {n}b_{n}(\mathbf{x}). \tag{23}\]  
These can be approximated by means of a numerical collocation scheme, i.e., by solving  
\[\int_{D}C(\mathbf{x}_{k},\mathbf{y})b_{n}(\mathbf{y})\mathrm{d}\mathbf{y}= \theta_{n}b_{n}(\mathbf{x}_{k}),\quad k=1,2,\ldots,M, \tag{24}\]  
in some well-chosen collocation points \(\mathbf{x}_{k}\). Following the [PERSON] method, see [33], the integral in Equation (24), is approximated by a numerical integration scheme which uses the collocation points as quadrature nodes, i.e.,  
\[\sum_{q=1}^{M}w_{q}C(\mathbf{x}_{k},\mathbf{y}_{q})\widetilde{b_{n}}(\mathbf{ y}_{q})=\widetilde{\theta}_{n}\widetilde{b_{n}}(\mathbf{x}_{k}),\quad k=1,2, \ldots,M. \tag{25}\]  
In matrix notation, this becomes  
\[\Sigma W\widetilde{B}_{n}=\widetilde{\theta}_{n}\widetilde{B}_{n}, \tag{26}\]  
where \(\Sigma\) is a symmetric positive semi-definite matrix with entries \(\Sigma_{k,q}=C(\mathbf{x}_{k},\mathbf{y}_{q})\), \(W\) is a diagonal matrix containing the weights \(w_{q}\) on its diagonal and \(\widetilde{B}_{n}\) is a vector with entries \(\widetilde{B}_{n,q}=\widetilde{b}_{n}(\mathbf{y}_{q})\). The matrix eigenvalue problem, Equation (26), can be reformulated as  
\[\Psi\widetilde{B}_{n}^{*}=\widetilde{\theta}_{n}\widetilde{B}_{n}^{*}, \tag{27}\]  
where \(\widetilde{B}_{n}^{*}=\sqrt{W}\widetilde{B}_{n}\) and \(\Psi=\sqrt{W}\Sigma\sqrt{W}\). Matrix \(\Psi\) is symmetric positive semi-definite. This implies that the eigenvalues \(\widetilde{\theta}_{n}\) are nonnegative real and the eigenvectors \(\widetilde{B}_{n}^{*}\) are orthogonal to each other. Using Equation (25), the Nystrom eigenfunctions \(\widetilde{b}_{n}(\mathbf{x})\) are obtained:\[\widetilde{b}_{n}(\mathbf{x})=\frac{1}{\widetilde{\theta}_{n}}\sum_{q=1}^{M}\sqrt{w _{q}}\widetilde{B}_{n,q}^{*}\mathcal{C}(\mathbf{x},\mathbf{y}_{q}), \tag{28}\]  
where \(\widetilde{B}_{n,q}^{*}\) stands for the q-th element of eigenvector \(\widetilde{B}_{n}^{*}\). These eigenvalues and eigenfunctions can be used as an approximate eigenpair in the KL expansion after a suitable normalization.