I think this post in the SAGA GIS forum might prove useful in answering your question about how slope is calculated:
https://sourceforge.net/p/saga-gis/discussion/354013/thread/27ecdc6b/
Also, based on my understanding of TWI (as a PhD hydrology student involved in hydrologic modeling), the D-Inf (Tarboton), MFD-md (Qin), DEMON (Costa-Cabral), and MFD (Quinn) with the exponent p=1.1 (Freemann) are the best options for determining the accumulation area 'a' in the TWI calculation.
I think Sorensen's work, and Qin's work, lends credence to my own semi-professional opinion. However, Qin's improved algorithm (MFD-md) hasn't been as throughly tested and used as much as the others have.
When I have used SAGA to calculate TWI, I first calculate the slope using the slope commend in the Terrain Analysis / Morphometry / Slope, Aspect, Curvature module, using the default algorithm mentioned in the forum post. Then I calculate catchment area using the Terrain Analysis / Hydrology / Catchment Area / Catchment Area (Parallel) choice, using either the MFD algorithm or the D-INF algorithm with 1.1 as the convergence factor (p=1.1, from Freeman).
Then I run the TWI option, in Terrain Analysis / Hydrology / Topographic Indices / Topographic Wetness Index (TWI), with the options to convert area to "1/cell size" and use standard calculation. I convert to specific catchment area because this is what the original formulation by Beven and Kirkby called for. As for the difference between the "Standard" and "TOPMODEL", I am not sure what it is - looking into that myself right now.
Reference page 106 onward of the linked pdf for more specific help: http://sourceforge.net/projects/saga-gis/files/SAGA%20-%20Documentation/SAGA%
20Documents/SagaManual.pdf/download
I forgot to add, this all assumes that the DEM has been preprocessed by filling sinks. That is another complicated topic, with two (main) separate options.
I hope this helps!
Tom
P.S. Edits for @reima:
This is something I have only recently dug into, and I can admit I don't think I've reached the bottom yet! I prefer the method of Lindsay and Creed, the minimum impact approach that chooses breeching or filling based on minimizing total topographical impact (officially named "Impact Reduction Algorithm" - IRA) which I thought was implemented in his Terrain Analysis Software tool (formerly TAS, now WhiteBox GAT - link: http://www.uoguelph.ca/~hydrogeo/Whitebox/).
However, even his tool seems to implement other filling schemes:
The sink/depression filling algorithm (basic, but insanely fast) by Wang and Liu (2006) - which I don't believe operates in the IRA manner, but similar to the way ArcMap fills sinks/depressions, just straight up without any breeching.
And the sink/depression filling of Planchon and Darboux (2001), which floods a DEM and then removes the water bit by bit - it can enforce a slope on the filed area, which I think might improve TI calculations.
ArcMap has a new "de-pitting" add-on (http://blogs.esri.com/esri/arcgis/2013/03/05/optimized-tool-for-dem-pit-removal-now-available/) that seems similar to Lindsay and Creeds IRA, but I haven't read the cited paper yet to determine how similar. This method might be worth a look.
I'm also interested in scrutinizing my assumption that TI calculations need filled DEMs. I have three different sized watershed DEMs (<100 sq km, 100-1000 sq km, >1000 sq km), clipped using a shape file from 10 m NED data. These are not filled, since the shape file already provided the watershed delineation. I am going to run the SAGA GIS TI calculation (MFD, p=1.1) on all three watersheds, on both filled and unfilled DEMs, using ArcMaps filling scheme (old and new), and the Wang and Liu algorithm (in Whitebox, maybe in SAGA), and the Planchon and Darboux algorithm (in Whitebox, maybe in SAGA). I will also be calculating the TI values using the TI calculation embedded in my hydrological model.
If you want, I can share these results with you. I might not have them for a month or so though, as I have other more pertinent research that my focus is currently on, but I need to refine my TI calculation process by mid May at the latest.
In spatial hydrology, DEM-based flow accumulation operations are typically static. That is, they represent a steady-state condition of the discharge of surface and near surface water passing through a point. Flow accumulation grids are actually accumulating contributing area downslope, i.e. areas not volumes. The accumulated cells that you described is because clearly the number of upslope pixels is related to the upslope area through the grid cell area. However, in addition to total number of upslope grid cells, it is common to represent this value as a 'specific contributing area' which is normalized for the contour length through the cell. Regardless of the units, the idea is that the contributing area can be used as a surrogate for the discharge, and for most locations this assumption is quite good in that there is a strong (often non-linear) relation between upslope area and discharge (at least for most streams).
However, you're talking about accumulating water volume and that requires a bit more information and likely implies that you would need to consider the dynamics of a full rainfall-runoff model. In particular, in addition to the upslope inputs (which the flow accumulation operation represent), you'd also need to account for the various losses in each grid cell due to evapotranspiration, percolation to deep groundwater, etc. This will be a function of much more than simply the topography, such as vegetation type and density, soil properties, the relation with the ground water table and much more. You may even need to account for spatial variation in the input of rainfall. Rainfall-runoff models usually require substantial parameterization to account for these factors. Thus, the topographically driven flow routing component represented by a flow accumulation grid would only be one such parameter (accounting for the redistribution of water within catchments) in a dynamic model which would then allow you to compute flow volumes. For an example of such a model, TOPMODEL is a good and easy-to-conceptualize starting point.
Best Answer
ArcMap uses the D8 routing algorithm for determining flow directions between cells, with no convergence exponent since it is a single direction algorithm, while the default in SAGA GIS (I believe - at least it is in mine) is MFD with the convergence you specify.
These methods differ and therefore will generate different catchment area results.
You could run SAGA GIS with the D8 algorithm (and a convergence = 1), and this should make the results match more closely, but I have the feeling that even then they will not be identical. Different implementations of the same algorithm (differences in the way they are coded up) should still lead to differences in the results.
Hope this helps.
Edit: See the comment from @WhiteBoxDev below as well.