Hydraulic fracturing processes are surrounded by uncertainty, as available data is typically scant. In this work, we present a sampling-based stochastic analysis of the hydraulic fracturing process by considering various system parameters to be random. Our analysis is based on the Perkins-Kern-Nordgren (PKN) model for hydraulic fracturing. This baseline model enables computation of high fidelity solutions, which avoids pollution of our stochastic results by inaccuracies in the deterministic solution procedure. In order to obtain the desired degree of accuracy of the computed solution, we supplement the employed time-dependent moving-mesh finite element method with two new enhancements: (i) global conservation of volume is enforced through a Lagrange multiplier; (ii) the weakly singular behavior of the solution at the fracture tip is resolved by supplementing the solution space with a tip enrichment function. This tip enrichment function enables the computation of the tip speed directly from its associated solution coefficient. A novel incremental-iterative solution procedure based on a backward-Euler time-integrator with sub-iterations is employed to solve the PKN model. Direct Monte-Carlo sampling is performed based on random variable and random field input parameters. The presented stochastic results quantify the dependence of the fracture evolution process-in particular the fracture length and fracture opening-on variations in the elastic properties and leak-off coefficient of the formation, and the height of the fracture. Keywords Hydraulic fracturing • Perkins-Kern-Nordgren model • Finite element method • Moving-boundary problem • Stochastic analysis • Random fields • Sensitivity analysis • Monte-Carlo method Hasini Garikapati