The problem of the precise number that lies on the surface is much more difficult because the number is quite random in a sense, and I don't know if much research has been done on it. If $a_r$ denotes the number of lattice points on the surface of the 3-d sphere with radius $r$ centered at the origin, then each individual $a_r$ fluctuates quite erratically.
Analysts need smoother behavior, so we study the sum instead: $ a_1 + a_2 + \cdots + a_r.$ This sum ends up being the number of lattice points contained inside and on the surface of the sphere of radius r, and hence this problem has been studied using analytic techniques successfully (see here ).
For your problem you require more algebraic flavored techniques. I'm not too sure about this type of work on this problem other than my own work here. As you'll soon find out reading the paper, my results are solely for the theoretical interest of having a generating function for the number of lattice points on the surface of a $N$-sphere, and practically using it is more computationally expensive than even brute force techniques. Hopefully you find some use or interest in it.