We implemented an extended version of the space-time ETAS model, which simultaneously incorporates earthquake focal depths and rupture geometries of large earthquakes and applied to the 2016 Kumamoto earthquake sequence. Results show that the new model corrects the estimation biases of model parameters due to the isotropic response function of aftershock locations adopted in the point source ETAS model, and that the reconstructed patterns of aftershock productivity density along the mainshock rupture plane show significant migrations in space and time along the mainshock fault. Another interesting phenomenon is that large aftershocks tend to nucleate at the edges of high productivity areas. The decay of direct aftershocks near the mainshock rupture is consistent with static stress changes caused by the mainshock. In simulations, the incorporation of focal depths improves the forecasting resolution. Through comparisons between aftershock productivity and slip distribution of the mainshock, we find that high aftershock productivity forms complementary patterns for coseismic slips in space, indicating that aftershocks play a role in the postseismic relaxation process.