The Triassic Shublik Formation of the Alaska North Slope is a world-class resource rock and has been identified as the major source of many of the conventional hydrocarbon accumulations on the North Slope, including Prudhoe Bay. Recent interest in the Shublik as a potential shale resource play has highlighted the need for robust hydraulic fracture modeling and simulation of the interval, but little geologic information is available because of the remote nature of the region and the complex character of the Shublik. In this study, a methodology was developed for identifying the critical variables needed for accurate planning of a hydraulic fracturing treatment in a play like the Shublik where much of the geology remains unconstrained. These identified critical variables can be used to develop a proxy model that can be used in lieu of a numerical simulator.This study was conducted in two stages. The first stage used 2-level fractional factorial design to identify the statistical significance of the input variables on the simulated fracture geometry. This stage was conducted in three phases, each phase incorporating progressively more complex assumptions about geology. Using the three most significant variables identified from first stage, the second stage of this study applies Box-Behnken experimental design and response surface methodology for quantifying functional relationships between input variables and the predicted fracture geometry. A pseudo 3D numerical simulator (Fracpro PT) and MATLAB were used to develop proxy models. These proxy models, typically a polynomial equation, are an easier alternative to Fracpro PT and can predict the fracture geometry with very less computational time.The use of experimental design drastically reduces the number of simulations required to evaluate large number of variables. With only 137 simulations, 26 variables were ranked based on their statistical significance and a non-linear proxy model was developed. Predicted values of the fracture geometry obtained using the proxy models were in good agreement with the simulated values of the fracture geometry (R 2 value of 99.39% for fracture length, R 2 value of 99.54% for fracture height and R 2 value of 98.17% for fracture width).