I've had some trouble with making a prototype (I don't seem to be much good with collision detection, even 2d collision is non trivial) so I thought I would outline the method and we can discuss whether it is worth pursuing.
@moopli I wanted to check what you are thinking about layers. So I was imagining the model would be 1 layer deep, basically there would be no movement in the radial direction, everything would be tangential. However you said "Other continental points would sit between and around cratons, in layers 1-2 points thick, to form full continents," which sounds rather different.
One of the issues about the spherical geometry is that points which start close together and which start moving parallel to each other can end up overlapping on the other side of the sphere. Now maybe that is not such a big deal (I imagine it would take a long time for a tectonic plate to orbit a sphere) but it is an issue, here is a diagram. So I'm not sure that simply interpolating quaternions will work very well, at least not globally.
So this is an algorithm which might work to avoid this problem (and allow for some nice dynamics).
0) Cover the sphere in a vector field (interestingly this runs into the hairy ball theorem and so will need to vanish somwehere). This represents the motion of the magma underneath the plates. It would also be possible to have a temperature parameter for the planets core, so it can cool over time and the magma can slow.
1) Plates are made of points and each point has a quaternion rotation axis and a speed, use this to move the points a small step.
2) For each point switch from a spherical geometry to a locally euclidean 2d patch around the point. That is take a gradient of the function which defines the sphere = (2x,2y,2z) and find it's value at the point of interest. Then work in the space with zero inner product with this vector.
3) Add up the forces in this tangent space, so if the point is alone far from any other it will just the the force of the magma pushing from below. However if the point is connected to others, or if it collides with others, forces will need to be added for these interactions to prevent the points from overlapping. This could be done by a usual collision detection system, I am having trouble implementing one.
Also there is the issue of solving the equations simultaneously, so if A collides with both B and C then we need to make sure if we run the algorithm for A then B then C that somehow this does not require the equations for A to be re run.
4) Convert the net force from this space back into an axis of rotation and go to step 1.
Anyway I am concerned this is pretty complicated. however it should provide some nice dynamics.