Hydraulic fracturing modelling is a multi-physics problem which must be solved in the reservoir domain, fracture domain, and wellbore domain. The physical phenomena that need to be quantitatively modelled include solid deformation and fluid flow in the poroelastic reservoir, fluid flow in the fractures, and fluid and proppant transport in the wellbore. These problems which are tightly coupled with each other in a highly non-linear manner are solved together implicitly for the first time. This fully implicit, parallelized 3D hydraulic fracturing simulator is capable of simulating simultaneous propagation of multiple fractures in multi-well pads. In this new simulator, four sets of governing equations are coupled implicitly and solved simultaneously. The solid deformation equations and reservoir pressure equations are discretized using the finite volume method while the fracture pressure equation is discretized using the finite area method. Fluid flow in the wellbore is modelled considering wellbore compressibility, perforation pressure drop, and wellbore friction pressure drop. A Barton-Bandis contact model is applied to model fracture closure due to stress shadow effects as well as production. The simulator is parallelized using the distributed memory approach with domain decomposition to improve the computation efficiency. We first validated our simulator with two well-known model problems and the numerical results show good agreement with the existing analytical solutions. We compare the performance of the new simulator with the past explicit methods and show that the new simulator offers significant improvement in numerical stability and computation speed. Finally, the simulator is applied to simulate and analyze one stage of a fracturing treatment with five fractures propagating simultaneously to show the unique capabilities of our simulator.