A two-dimensional simulation modeling that has been performed in a self-consistent way for analysis on the fully coupled transports of plasma, recycling neutrals, and intrinsic carbon impurities in the divertor domain of tokamaks is presented. The numerical model coupling the three major species transports in the tokamak edge is based on a fluid-particle hybrid approach where the plasma is described as a single magnetohydrodynamic fluid while the neutrals and impurities are treated as kinetic particles using the Monte Carlo technique. This simulation code is applied to the KSTAR ͑Korea Superconducting Tokamak Advanced Research͒ tokamak ͓G. S. Lee, J. Kim, S. M. Hwang et al., Nucl. Fusion 40, 575 ͑2000͔͒ to calculate the peak heat flux on the divertor plate and to explore the divertor plasma behavior depending on the upstream conditions in its base line operation mode for various values of input heating power and separatrix plasma density. The numerical modeling for the KSTAR tokamak shows that its full-powered operation is subject to the peak heat loads on the divertor plate exceeding an engineering limit, and reveals that the recycling zone is formed in front of the divertor by increasing plasma density and by reducing power flow into the scrape-off layer. Compared with other researchers' work, the present hybrid simulation more rigorously reproduces severe electron pressure losses along field lines by the presence of recycling zone accounting for the transitions between the sheath limited and the detached divertor regimes. The substantial profile changes in carbon impurity population and ionic composition also represent the key features of this divertor regime transition.