This study proposes a modeling method for heterogeneous rock masses using the discrete element method (DEM) and investigates the seepage characteristics of heterogeneous rock masses. The Weibull distribution is used to describe the heterogeneity of rock hydraulic properties and integrated into the Fish program, which is based on the DEM. This program overcomes the limitation of the Universal Distinct Element Code (UDEC) software, which cannot exceed 50 parameter groups. A method for parameter assignment of heterogeneous rocks is proposed, along with a method for estimating the initial flow rate value of heterogeneous models. Based on the established heterogeneity calculation model, the influence of block homogeneity, hydraulic aperture homogeneity, and stress on the seepage characteristics is studied. The results indicate that under zero stress conditions, the flow rate is positively correlated with N^0.5, showing a strong linear relationship. The linear relationship is gradually enhanced with the increase in the shape parameters. The relationship between the flow rate and shape parameters is logarithmic with a correlation coefficient greater than 0.9654. The relationship between the flow rate and the axial pressure and confining pressure can be described by quadratic and cubic polynomials, respectively. The study further discusses the variation characteristics of equivalent hydraulic apertures under various axial pressures, confining pressures, and shape parameters. The results show that the total flow is positively correlated with N^0.5 and shows a strong linear relationship. The relationship between the total flow and the shape parameter is logarithmic with a correlation coefficient greater than 0.9654. The relationship between the flow rate and the axial pressure is a quadratic polynomial with a correlation of more than 99.72%. The relationship between the flow rate and the confining pressure is a cubic polynomial with a correlation of more than 99.98%. In identical stress environments, the increase of the shape parameters corresponds to a gradual increase in the seepage flow rate. The study also shows that the equivalent hydraulic aperture (EHA) linearly depends on the confining pressure and has a quadratic polynomial relationship with the axial pressure. The EHA has a logarithmic relationship with the shape parameters. The study concludes that the proposed method is reasonable, feasible, and effective for modeling heterogeneous rock masses. The results provide a theoretical basis for understanding the seepage characteristics of heterogeneous rock masses and have practical significance for engineering applications.This study proposes a modeling method for heterogeneous rock masses using the discrete element method (DEM) and investigates the seepage characteristics of heterogeneous rock masses. The Weibull distribution is used to describe the heterogeneity of rock hydraulic properties and integrated into the Fish program, which is based on the DEM. This program overcomes the limitation of the Universal Distinct Element Code (UDEC) software, which cannot exceed 50 parameter groups. A method for parameter assignment of heterogeneous rocks is proposed, along with a method for estimating the initial flow rate value of heterogeneous models. Based on the established heterogeneity calculation model, the influence of block homogeneity, hydraulic aperture homogeneity, and stress on the seepage characteristics is studied. The results indicate that under zero stress conditions, the flow rate is positively correlated with N^0.5, showing a strong linear relationship. The linear relationship is gradually enhanced with the increase in the shape parameters. The relationship between the flow rate and shape parameters is logarithmic with a correlation coefficient greater than 0.9654. The relationship between the flow rate and the axial pressure and confining pressure can be described by quadratic and cubic polynomials, respectively. The study further discusses the variation characteristics of equivalent hydraulic apertures under various axial pressures, confining pressures, and shape parameters. The results show that the total flow is positively correlated with N^0.5 and shows a strong linear relationship. The relationship between the total flow and the shape parameter is logarithmic with a correlation coefficient greater than 0.9654. The relationship between the flow rate and the axial pressure is a quadratic polynomial with a correlation of more than 99.72%. The relationship between the flow rate and the confining pressure is a cubic polynomial with a correlation of more than 99.98%. In identical stress environments, the increase of the shape parameters corresponds to a gradual increase in the seepage flow rate. The study also shows that the equivalent hydraulic aperture (EHA) linearly depends on the confining pressure and has a quadratic polynomial relationship with the axial pressure. The EHA has a logarithmic relationship with the shape parameters. The study concludes that the proposed method is reasonable, feasible, and effective for modeling heterogeneous rock masses. The results provide a theoretical basis for understanding the seepage characteristics of heterogeneous rock masses and have practical significance for engineering applications.