The maximum possible size is the number of mesh volume elements times the number of geometry cells (the track lengths do not need to be stored, as I was reminded at the CNERG meeting). You could allocate this all ahead of time using record arrays but this seems wasteful. Record arrays would also allow volume fractions to be stored side-by-side with their standard errors. 

I like scipy.sparse arrays because it seems like they would have smaller memory footprints. The a minor disadvantage is in order to store the results and the standard errors, either a 3D array or two 2D arrays are needed. If two 2D arrays are used, the results and standard errors could be returned separately, which wouldn't be a bad idea considering a lot of use cases may not utilize the standard errors.

Both seem easy enough to implement, though I think we should decide now. I think an optimized sparse array makes more sense then an non-optimized sparse record array, or a record array that changes size.