This research discusses the hydraulic fracturing process in a heterogeneous medium. To model this process, the PPR- potential-based cohesive zone model (Park et. al.,2009; Park et. al., 2012) was implemented in the finite element method (FEM) to model the crack propagation; and the lattice-Bolztmann model (LBM) was used to represent the fluid injected into the fracture. These models are coupled in an iterative process (Mejía Camones et al., 2014; Mejía Camones et. al., 2015). In the FEM, the crack’s boundary is represented by the faces of cohesive elements. The position and velocity of these faces are transferred to the LBM as a boundary condition, where the forces applied on these boundaries, caused by the fluid injected, can be calculated. These forces are transferred to the FEM as external force applied on the element's faces, where the new position and velocity of the boundary of the crack is computed. These values are transferred to the LBM, initializing a new calculation cycle. The utilization of cohesive elements allows introducing different energy and strength fracture values in each element. Thus, the crack propagation in a heterogeneous medium can be modeled. This research shows the results of the crack propagation in heterogeneous medium, comparing with experimental results.