Choosing the right algorithm for hydrologically correcting a DEM really depends on your particular application. DEM size is certainly an issue, as you discovered. If you have a massive DEM, there are only a few algorithms that will work for your application. Another important consideration is whether or not all of the topographic depressions in your DEM are actually artifacts and whether any of the real depression that may be in your DEM are significant in terms of the process that you are modelling (modelling surface flow in your case). If there are real features, and they are significant, you should look for a way to retain these features in your DEM. Then there is the depression filling versus depression breaching question. It turns out that although there are many depression filling algorithm available, they will all result in the same solution (see this question). Importantly depression filling almost always results in far greater impact to the DEM (i.e. more modification of the original elevation values) than does depression breaching. I almost always recommend depression breaching solutions over filling.
So now for the really tricky part. You asked for a tool to recommend. As full disclosure, I am the lead developer of the open-source GIS Whitebox Geospatial Analysis Tools. Whitebox GAT has several tools available for depression removal within the DEM Pre-processing toolbox contained within the Hydrological Tools toolbox. For example, there are tools that implement the Planchon and Darboux algorithm, which is suitable for use with massive DEMs, as well as the Wang and Liu (2006) filling algorithm, which is far faster but not suitable for massive DEMs. It also contains tools for depression breaching, which is a very good solution for fine-resolution massive DEMs. I've just finished writing a tool that performs depression breaching very efficiently (as fast as the Wang and Liu algorithm) but it is not as suited to massive DEMs. So there is a lot of options there. Welcome to the murky waters of DEM pre-processing! Many people, I think, just use whatever tool is available in ArcGIS without much thinking about how greatly it can impact the results of the analysis.
EDIT:
Okay, if it's channel lines that you're after then there are really two general approaches to accomplishing this task. The first approach is to use a morphologically based approach. These techniques are based on the early work of Peucker and Douglas (1975) and Johnston and Rosenfeld (1975), and are based on finding v-shaped profiles, some measure of thresholded curvature, or some low topographic position. They have the advantage that they don't require depression filling or removal of flat areas and they can work quite well for highly accurate, fine resolution DEMs in areas of relatively high topographic relief. They have the disadvantage that they don't necessarily produce continuous, single-cell wide channel networks and therefore require substantial post-processing. They can also be greatly affected by the scale of operation as seen by the impacts of DEM smoothing filters. In Whitebox GAT, you will find the Find Valleys series of tools in the Stream Network Analysis toolbox that can be used to perform this type of channel mapping.
The second general approach is based on measuring upslope contributing area. These techniques all go back to the work of O'Callaghan and Mark and the D8 flow algorithm, although there have been numerous studies since. The general idea is that you can use the DEM as a means of discerning flowpaths and for measuring upslope contributing area. In reality the importance of upslope contributing area is that it is a surrogate for discharge (the larger the area draining to a site, the greater the amount water discharged from the location). While drainage area and discharge are actually quite well correlated for most basins, it is generally a non-linear relation. With the DEM handling both the flow direction problem and the upslope area/discharge measure, the issue becomes one of defining channel heads, i.e. it is a channel initiation issue related to erosion of surface materials. We generally figure this out simply by applying a constant landscape invariant threshold to the contributing area map, although there is a fair amount of research that demonstrates that no such constant threshold actually exists and that it is dependent on local slope and surface materials (see the work of Montgomery). That's why you generally have to experiment a fair amount with various threshold before finding one that's appropriate. Anyhow, in Whitebox GAT, you can use one of the several flow routing algorithms available in the Hydrological Analysis toolbox to measure upslope area, although remember that these techniques DO require removal of depressions and flat areas to work properly. Also see the Extract Streams tool in the Stream Network Analysis toolbox. The main advantage of these area based methods is that they can produce continuous, single-cell wide channels networks that are amenable to stream network analysis. The tools in the Stream Network Analysis toolbox can be used to do all sorts of analytical processing of raster stream networks extracted from DEMs including stream classification and ordering, channel slope analysis, profile mapping, identification of main channels, etc. Importantly for you, there is a tool called Raster Streams to Vector that will convert your raster stream into a shapefile that you can then use for display purposes.
