A conformal projection on a world scale is a good idea, but it won't get the slopes right. Conformality assures that the scale distortion at any point does not vary with direction: that makes the slopes in all directions comparable, resulting in relatively accurate computation of aspect and of flow drainage directions. However, all conformal maps covering large areas must introduce substantial distortion from one point to the other, which means they cannot be used by themselves to compute a consistent slope across the entire study area.
One solution is to correct the computed slopes. This method is described in my answer at Global DEM to slope calculation.
Having done this, perform the followup calculation of flow directions, aspects, etc. using the corrected slopes. Then, for carrying out a flow accumulation calculation (or other calculations that depend on slopes and aspects), consider reprojecting those DEMs to an equal-area projection. There are many global equal-area projections to choose from; for instance, any cylindrical equal area projection will be fine for latitudes up to 60 degrees (or even almost to the poles).
One potential problem to anticipate concerns the resampling of values that occurs during reprojection. Most flow accumulation algorithms discretize the flow directions. For instance, ArcGIS uses a "D8" algorithm in which only eight directions are recorded, employing a binary code comprising powers of two from 1 = 2^0 through 128 = 2^7. Terrible things can happen to these values if they are averaged during reprojection. For instance, direction 127 is next to direction 1: during resampling just about any direction between 1 and 127 can result near those locations, which can create a huge mess. One also has to wonder what artifacts a nearest-neighbor resampling procedure (which doesn't average) would produce.
I have not experimented with the following solution and so can only propose it as an idea for further consideration and testing. As with most "circular" data, an effective approach to averaging is to replace directions by their cosines and sines (essentially converting a polar coordinate into Cartesian coordinates). This could be done for flow direction codes merely with a lookup table. For example, ESRI software codes would be transformed as follows:
Code Direction Bearing Cosine Sine
---- --------- ------- ------ ----
1 E 90 0 1
2 SE 135 -s s
4 S 180 -1 0
8 SW 225 -s -s
16 W 270 0 -1
32 NW 315 s -s
64 N 0 1 0
128 NE 45 s s
where s = Sqrt(1/2) = 0.7071... . (I would be tempted to multiply these values by the tangent of the slope, turning them into components of a vector field representing the gradient of the surface. Reprojection would then be averaging the gradient itself.)
After reprojection, the resulting (cosine, sine) pairs can be converted back into flow direction codes according to the same scheme. The inverse tangent of the interpolated (cosine, sine) determines an angle which can be written in degrees between 0 and 360. Dividing by 45 = (360/8) gives a number between 0 and 8. Rounding this to the nearest whole number determines the flow direction. For example, suppose the interpolated cosine equals -0.1 and the interpolated sine equals 0.9. (Notice their squares do not sum to unity! The separate interpolations will usually destroy this relationship.) The inverse tangent of this pair is 96 degrees which, when divided by 45, equals 2.14. That rounds to 2, corresponding to 2 * 45 = 90 degrees, for a direction code of 1.
The same approach will work with any flow direction coding scheme: convert the codes into (cosine, sine) pairs, reproject, then convert back after determining the angle via the inverse tangent.
For an assessment of the distortions in flow direction that arise when these issues are not addressed, please read the reply at Calculating flow direction and delineating basins....