This paper describes a new fast and scalable parallel algorithm to compute global flow accumulation for automatic drainage network extraction in large digital elevation models (DEM for short). Our method uses the D8 model to compute the flow directions for all pixels in the DEM (except NODATA and oceans). A parallel spanning tree algorithm is proposed to compute hierarchical catchment basins to model the flow of water from a sink (local minima) moving on DEM to its outlet (ocean, NODATA, or border of DEM). And finally, based on local flow accumulation and the hierarchical trees between sinks, we determinate entirely the global flow accumulation. From that, the drainage networks of DEM can be extracted. Our method does not need any preprocessing like stream burning on the initial DEM and tends to make the most of incomplete DEMs. Our algorithms are entirely parallel. Efficiency and scalability have been tested on different large DEMs.