## 1 Introduction  
There is an increasing need to accurately simulate and compute solutions to engineering problems, whilst accounting for model uncertainties. Methods for uncertainty quantification and propagation in engineering disciplines can be categorized into two groups: non-sampling and sampling methods. Examples of non-sampling methods are the perturbation method and the Stochastic Galerkin Finite Element method. The perturbation method is based on a Taylor series expansion approximating the mean and variance of the solution [1]. While quite effective, its use is restricted to models with a limited number of relatively small uncertainties. The Stochastic Galerkin method, proposed by [PERSON] and [PERSON] [2], is based on a spectral representation in the stochastic space. It transforms the uncertain coefficient partial differential equation (PDE) problem by means of a Galerkin projection technique into a large coupled system of deterministic PDEs. This method allows for somewhat larger numbersof uncertainties and is quite accurate. However, it is highly intrusive and memory demanding, making its implementation cumbersome, and restricting its use to still rather low stochastic dimensions. Sampling methods, on the other hand, are typically non-intrusive. Each sample corresponds to a deterministic solve, with a particular choice for the parameter values. Two popular examples are the Stochastic Collocation method, see [3], and the Monte Carlo (MC) method, see [4]. The former samples a stochastic PDE at a carefully selected multidimensional set of collocation points. After this sampling, a Lagrange interpolation is performed, leading to a polynomial response surface. From this response surface, the relevant stochastic characteristics can easily be computed in a post-processing step. Like Stochastic Galerkin, the Stochastic Collocation method suffers from the curse of dimensionality: the computational cost grows exponentially with the number of random variables considered in the problem. The MC method on the other hand, selects its samples randomly and does not suffer from the curse of dimensionality. However, a drawback is its slow convergence as a function of the number of samples taken. The convergence of Monte Carlo can be accelerated in a variety of ways. For example, an alternative non-random selection of sampling points can be used, as in Quasi-Monte Carlo [5, 6] and Latin Hypercube [7] sampling methods. Other techniques, such as Multilevel Monte Carlo (MLMC) [8], Multilevel Quasi-Monte Carlo (MLQMC) [9] and its generalizations, see, e.g., [10, 11], can speed up the method. These improved Monte Carlo methods are based on a hierarchy of meshes with increasing resolution, where samples on coarser meshes are computationally less expensive than samples on finer meshes. For completeness, we mention that there also exist hybrid variants which exhibit both a sampling and non-sampling character. This type of methods combine, for example, the Stochastic Finite Element methodology with Monte Carlo sampling or a multi-dimensional cubature method, see, e.g., [12, 13]. In previous work, we presented a comparison between Monte Carlo and Multilevel Monte Carlo [14], and between Multilevel Monte Carlo and Multilevel Quasi-Monte Carlo [15]. In this work, we will present a novel multilevel method: p-refined Multilevel Quasi-Monte Carlo (p-MLQMC). This method reduces the cost of the simulation by taking samples on a hierarchy of nested p-refined meshes, and by selecting the sample points in a non-random way. In [16], the authors introduced a multi-order Monte Carlo algorithm (MOMC) for computing the statistics of quantities of interest described by linear hyperbolic equations with stochastic parameters. The MOMC method relies on a hierarchy of p-refined meshes, similar to the setup in this paper. A numerical analysis of MOMC, including numerical results for dispersion in long-distance propagation of waves subjected to uncertainty can be found in [16].