Hydraulic fracturing has gained increasing attention as it allows the constrained natural gas and crude oil to flow out of low-permeability shale formations and significantly increase production. Perilous operating states of extremely high pressure also raise some safety concerns, requiring us to formulate an appropriate dynamic model, and provide a careful engineering control to ensure safe operating conditions. Moreover, uncertainties due to spatially varying rock properties increase the difficulties in control of the fracturing process. In this work, we formulate a first-principles model by considering the fracture evolution, mass transport of substances in the slurry, changing fluid properties, and the monitored operating pressure on the ground level.Next, we implement nonlinear model predictive control (NMPC) to control the process under a set of final requirements and process constraints. Our results show that the performance of standard NMPC degrades when the rock uncertainty causes the parameter mismatch between the process and the predictive model in the controller.With standard NMPC, designed with a nominal model, the process fails to meet the terminal requirements of fracture geometry, and pressure is violated in one of the parameter mismatch cases. Therefore, we resort to multistage NMPC, which considers uncertainty evolution in a scenario tree with separate control sequences to address constraint violations. We demonstrate that multistage NMPC presents good performance by showing constraint satisfaction whether the uncertain rock parameter realization is time-invariant or time-variant. We also simulate the process with multistage NMPC including different numbers of scenarios and compare their control performance. Our investigation demonstrates that multistage NMPC effectively manages parametric uncertainties attributed to non-homogeneous rock formation, and provides a promising control strategy for the hydraulic fracturing process.