The vertical variation in elastic moduli in layered rocks has been neglected when the displacement discontinuity method (DDM) is applied to hydraulic fracturing simulations. This is primarily due to theoretical limitations for obtaining solutions to such problems. A new method for the computation of kernels in layered media based on the method of images is presented in this paper, which allows us to apply the DDM to layered rock formations. It can evaluate the contribution of multiple images in contrast to previous researches using bi-material solutions. The new method was verified by comparing kernels with analytical solutions. It is applied to crack problems in three-layered media, and the effect of element types, the number of images, and the number of Gaussian points are investigated. It is found that bi-material solutions are sufficient for the calculation of the width of a single crack in three-layered media under the plain strain condition.