Excessive water production in mature heavy oil fields causes incremental costs, energy consumption, and inefficiency. Understanding multiphase flows near the wellbore is an alternative to improve production efficiency. Therefore, this study conducts a series of numerical experiments based on the full set of the Navier-Stokes equations in 3D to simulate multiphase flows in porous media for heavy oil production horizontal wells. The solution given by this advanced mathematical formulation led to the description of the movement of the fluids near the wellbore with unprecedented detail. A sensitivity analysis was conducted on different rock and fluid properties such as permeability and oil viscosity, assuming homogeneous porous media. The influence of these parameters on the prediction of the breakthrough time, aquifer movement, and the severity of water production was noticed. Finally, the numerical model was verified against field data using two approaches. The first one was conducting a history match assuming homogeneous rock properties. In contrast, the second one used heterogeneous rock properties measured from well logging, achieving a lower deviation than field data, about 20%. The homogeneous numerical experiments showed that the breakthrough occurs at the heel with a subsequent crestation along the horizontal well. Moreover, at adverse mobility ratios, excessive water production tends to happen in water connings at the heel with an inflow area less than 1% of the total inflow area of the completion liner. Different aquifer movement dynamics were found for the heterogeneous case, like the breakthrough through multiple locations along the horizontal well. Finally, critical hydraulic data in the well, such as the pressure and velocity profiles, were obtained, which could be used to improve production efficiency. The numerical model presented in this study is proposed as an alternative to conducting subsurface modeling and well designs.