Simulating the surface flow dynamics is of great importance in disaster prevention and mitigation, which could benefit the land remediation, regional planning, environmental protection and water resource management. Therefore, in this paper the Compute-Unified-Device-Architecture (CUDA) is embedded into the TIN-based surface flow dynamic simulation algorithm to get a parallel method for simulating the surface flow discharge. The aim of this study is to rapidly and accurately perform the surface flow dynamic simulation at any position and any time to meet the actual application requirements. First, the proposed algorithm extracts the critical points and drainage network from a high-precision Digital Elevation Model (DEM) to generate a drainage-constrained Triangulated Irregular Network (TIN). Second, the flow direction of each triangular facet over the TIN is calculated by the coordinate data of the triangular facets. Third, the flow path network is traced by the flow direction and rainfall source points. Fourth, the flow velocity of rain dropsover the flow paths is obtained by the flow velocity calculation kernel function based on the manning formula. Finally, by combing the flow velocity with the pre-set time, the algorithm can rapidly simulate the flow discharge at any location by usingthe flow discharge count kernel function. The specific paralle lization process consists of the transmission mode of data, partition model of thread and implement ation of the flow velocity calculation kernel function and flow discharge count kernel function. Data transmission in the paralle lization process includes the input and output of data. It inputs the data of the DEM, rainfall source points and flow path network from the CPU to GPU and outputs the data of the flow discharge calculated by the above two kernel functions from the GPU to CPU. The two kernel functions are parallelized by the flow paths. Each thread handles a single flow path. As a result, flow paths are divided by the partitioning method to obtain numerous thread blocks and the number of the thread in each thread block is allocated by the computing power of the GPU. The Black Brook Watershed (BBW) located in the north-eastern of New Brun swick was selected as the study area. To validate its accuracy and efficiency, the proposed method was compared with TIN-based surface flow dynamic simulation method. The experimental results demonstrate that the proposed algorithm can greatly improve the simulation efficiency while maintaining the accuracy simultaneously and its acceleration ratio can reach up to 11.2. In addition, the parallel algorithm was compared with the Soil and Water Assessment Tool (SWAT) model to verify its precision. The experimental results indicate that the Nash coefficient of the parallel method is increased by 88% and the correlation coefficient is increased by 5%.