Hydraulic fracturing has become a widely used method for oil/gas extraction across the whole world, especially for unconventional formations. However, the study of hydraulic fracturing is still limited to simple fracture networks. The existing models rarely consider fracture propagation and interactions of multiple fractures in complex fracture networks. In this thesis, the extended finite element method (XFEM) will be used to model the evolution of hydraulic fracturing occurring in specific conditions. In order to take the advantage of XFEM and make a further contribution to complex hydraulic fracturing, an effective solid-fluid coupling method is proposed to improve the coupling effect. A new scheme for multiple hydraulic fracturing is developed to study the hydraulic fracturing in complex fracture networks. XFEM offers great advantage in the study of fracture propagation as it avoids remeshing and mesh refinement, thus providing considerable saving in computational effort. In addition, XFEM can simulate crack propagation without pre-setting the path. It is shown that the commonly used narrow band level set method has drawbacks in tracking the crack path, resulting in inaccurate tracking outcomes. To address this shortcoming, a geometric tracking technique is proposed to efficiently track the crack path and search for enriched elements and nodes in 2D models. The study of multiple cracks and crack intersection in XFEM has not been well developed. Many software packages such as ABAQUS have strong limitations in modelling crack intersection and poor convergence for multiple cracks. In this research, based on the study of a single crack, the formulation of multiple cracks is presented by considering incorporation of enrichments of different cracks into one matrix to solve global stiffness matrix. In addition, the issue of crack intersection is analysed based on the intersecting geometry of two cracks and the redistribution of Gauss points. Two methods are proposed to improve the accuracy of computation of symmetric nodal displacements along cracks resulting from the approximation errors in XFEM. Furthermore, a method is presented to determine the stress intensity factor near the tip for the cracks subjected to surface loads. Several examples are used to explain the application of proposed methods. The modelling results are verified by comparing to analytical solutions and experimental tests. The research in this thesis will provide useful information for further study of complex hydraulic fracturing using XFEM.