The notion that the center of mass should correspond to "not rolling"
implies that the answer should be consistent with an embedding of
$S^2$ in $\mathbb R^3$.
(After all, if $S^2$ is the only space that exists, how could
the sphere "roll" anywhere?)
In fact, we be able to define the center of mass using that embedding,
and in particular it seems that this definition must place the
center of mass on the line that passes through the center of the sphere
and through the $\mathbb R^3$-center of mass of the points
(the center of mass of the corresponding collection of points $\mathbb R^3$).

The tricky part is that if we say that a collection of
points in $S^2$ with masses somehow acts like a single point with mass
$\sum m_i$ at point $\overline p$ in $S^2$, for example for the purpose
of finding the new center of mass when this set of points is combined
with another set (which works fine for the $\mathbb R^3$-center of mass),
we may get inconsistent results;
for example, given $n$ points with masses,
if we select $n-1$ points, replace them with a point of their
combined mass at their center of mass,
and then take the center of mass of that point and the $n$th point
from the original set, the result may depend on which $n-1$ points we
select first.
The reason is that projecting the $\mathbb R^3$-center of mass
onto the surface of the sphere increases the moment arm of the
"center-of-mass" point by some amount that depends on how much the $n-1$ selected points cancel each others' moment arms in the first place.
In other words, if center-of-mass calculations are considered as
weighted averages, the effect of the projection can "overweight"
a subset of the points when the points in that subset are combined together.

So my proposal is this:
embed $S^2$ in $\mathbb R^3$ so that the center of $S^2$
is at the origin of $\mathbb R^3$.
Compute $\overline p_{\mathbb R^3}$,
the $\mathbb R^3$-center of mass of the $n$ points.
If $\| \overline p_{\mathbb R^3} \| = 0$,
let $\overline p$ be any arbitrary point on $S^2$;
otherwise let
$\overline p = \frac{\overline p_{\mathbb R^3}}
{\| \overline p_{\mathbb R^3} \|}$.
Assign the mass
$\| \overline p_{\mathbb R^3} \| \sum m_i$ to $\overline p$.

The factor $\| \overline p_{\mathbb R^3} \|$ in the formula for
combining the masses of a subset of some points can be considered a
correction of the possible "overweighting" of that subset.

In $\mathbb R^3$, the moment of the mass
$\| \overline p_{\mathbb R^3} \| \sum m_i$ at $\overline p$
around any axis through the origin is equal to the moment around the
same axis of the mass $\sum m_i$ at $\overline p_{\mathbb R^3}$,
so when we combine a set of points with masses into a new point on $S^2$
with a new mass in this manner, the resulting point has the same effect on
new calculations of the center of mass as the original set of points,
just as it does when we compute centers of mass in $\mathbb R^3$.