Thanks for the replies. It does seem tricky and I'll carry on searching.
I started my model from cold with the rocky core at 100K and the
surrounding ice at 0K. Clearly wrong but a good starting point for
working with the model. As the planet heats up the T^4 radiation from
the surface seems to pin the surface to the same temperature as
surrounding space. Once my ocean forms I need a formula for the power
fow. For conduction I have:
def conductive_power_flow(temperature_1, temperature_2, area, thickness,
conductivity):
return (temperature_1-temperature_2) * area * conductivity / thickness
An example of a render of the system is attached (done with a free open
source package called Povray). I thought I'd create a movie too.
Regards,
Roy