Sophisticated catchment runoff problems necessitate conjunctive modeling of overland flow and subsurface flow. In this paper, finite difference numerical methods are studied for simulation of catchment runoff of two-dimensional surface flow interacting with three-dimensional unsaturated and saturated sub-surface flows. The equations representing the flows are mathematically classified as a type of heat diffusion equation. Therefore, two-and three-dimensional numerical methods for heat diffusion equations were investigated for applications to the surface and sub-surface flow sub-models in terms of accuracy, stability, and calculation time. The methods are the purely explicit method, Saul'yev's methods, the alternating direction explicit (ADE) methods, and the alternating direction implicit (ADI) methods. The methods are first examined on surface and sub-surface flows separately; subsequently, 12 selected combinations of methods were investigated for modeling the conjunctive flows. Saul'yev's downstream (S-d) method was found to be the preferred method for two-dimensional surface flow modeling, whereas the ADE method of Barakat and Clark is a less accurate, stable alternative. For the three-dimensional sub-surface flow model, the ADE method of Larkin (ADE-L) and Brian's ADI method are unconditionally stable and more accurate than the other methods. The calculations of the conjunctive models utilizing the S-d surface flow sub-model give excellent results and confirm the expectation that the errors of the surface and sub-surface sub-models interact. The surface sub-model dominates the accuracy and stability of the conjunctive model, whereas the sub-surface sub-model dominates the calculation time, suggesting the desirability of using a smaller time increment for the surface sub-model.