Fractals don't have surfaces and it's imposible to use standard 3D rendering techniques.
Just to point out that "Fractals don't have surfaces" is not strictly correct - they do have surfaces it's just that the true (infinite limit) surface is purely theoretical.
When rendering 3D/4D fractals the render always compromises by limiting the max. iteration depth and a given fractal formula at a given iteration depth does have a concrete, finite surface defined by |In(f(z))|=bailout i.e. the surface is where the magnitude of the nth iterate of f(z) is equal to the bailout value - the whole solid being defined by |In(f(z))|<=bailout.
As to ways of calculating the normals, I simply use the "found" surface adjacent "found" pixels method (or better still half-pixel offsets) that's good but doesn't show any sub-pixel detail - for that it's better to use a method such as the one you employ or to use the iteration density at pixel or sub-pixel offsets (or to find "solid" based on DE instead of iteration and do the same thing with the DE at pixel/sub-pixel offsets).
I'm still attempting to find a killer algorithm that will allow solid based on any colouring algorithm, see:
If no image above:http://makinmagic.deviantart.com/art/Alien-Interface-ANIM-29854835
That's a Julibrot solid on orbit trapping but in order to get (almost) visually flawless results I had to severely limit the max. iteration count (<10), I want to develop an algorithm that will render something similar but at say 60 iterations or even higher visually "accurately" and in reasonable times but without having to use masses of memory.
These attempts at pure raycasting orbit-trap heightfields illustrate the problem:http://makinmagic.deviantart.com/art/Orbit-Trap-Heightfield-128850151http://makinmagic.deviantart.com/art/Orbit-Trap-Height-Detail-128912127