The numerical modeling of hydraulic fractures in unconventional reservoirs presents significant challenges for field applications. There remains a need for accurate models that field personnel can use, yet remains consistent to the underlying physics of the problem [1]. For numerical simulations, several authors have considered a number of issues: the coupling between fracture mechanics and fluid dynamics in the fracture [2], fracture interaction [3][4][5], proppant transport [6], and others [7][8][9]. However, the available literature within the oil and gas industry often ignores the importance of the crack tip in modeling applications developed for engineering design. The importance of accurate modeling of the stress induced near the crack tip is likely critical in complex geological reservoirs where multiple propagating crack tips are interacting with natural fractures. This study investigates the influence of various boundary element numerical techniques on the accuracy of the calculated stress intensity factor near the crack tip and on the fracture profile, in general. The work described here is a part of a long-term project in the development of more accurate and efficient numerical simulations for field engineering applications.For this investigation, the authors used the displacement discontinuity method (DDM). The numerical technique is applied using constant and higher-order elements. Further, the authors also applied special crack tip elements, derived elsewhere [10], to capture the square root displacement variation at the crack tip, expected from Linear Elastic Fracture Mechanics (LEFM). The authors expect that special crack tip elements will provide the necessary flexibility to choose other tip profiles. The crack tip elements may prove instrumental for efficient modeling of the different near-tip displacement profiles exhibited by Viscosity-Dominated or Toughness-Dominated regimes in hydraulic fracture propagation. As others have shown [1,4,7], the accuracy of tip asymptote is critical in characterizing the stresses in the near-tip region of a propagating fracture.The authors examined the numerically derived stress intensity factor for several crack geometries with and without higher-order elements and with and without special tip elements, to analytical solutions. As expected, they found that the cases with higher-order elements and special tip elements provide more accurate results than the cases with constant displacement discontinuity and/or no tip elements. However, the numerical technique developed still proved efficient.These results show that numerical simulators can incorporate proper crack-tip treatments effectively. In addition, higher-order elements increase computational efficiency by reducing the number of elements instead of simply increasing the discretization of constant displacement elements. The accurate modeling of stress intensity factors is necessary to better simulate curved fractures, kinked cracks and interaction between fractures.