Numerous numerical models have been developed to simulate the hydraulic fracturing process in the subsurface, but few of them can predict hydraulic fracture branching, a key phenomenon producing complex fracture networks that enable economic hydrocarbon production from unconventional reservoirs. In this study, we present a 3D numerical model that combines fracture mechanics and poroelasticity to predict fracture complexity during hydraulic fracturing. In our model, hydraulic fracture propagation is governed by the fracture-mechanics-based Microplane 7 model and flow in fractures is governed by the cubic law with additional consideration of fluid leak-off at fracture walls. To solve the proposed model, we couple an in-house explicit finite element code, HOSS, for stress calculations and an open-source implicit finite volume program, PFLOTRAN, for fluid flow solutions. We validate the coupled model using Sneddon's analytical solution for fracture aperture distribution. Through three modeling cases, we discussed the hydraulic fracture branching mechanism: fluid leak-off at hydraulic fracture walls and the changes of Biot effective stress coefficient in the transverse direction after the primary hydraulic fracture crosses a pre-existing natural fracture provides a tensile body force for hydraulic fracture branching. The proposed model enables us to identify and optimize the key parameters that govern hydraulic fracture branching in a three-dimensional space.