A terrestrial planet is molten during formation and may remain so if subject to intense insolation or tidal forces. Observations continue to favour the detection and characterisation of hot planets, potentially with large outgassed atmospheres. We aim to determine the radius of hot Earth-like planets with large outgassed atmospheres and explore differences between molten and solid silicate planets and their influence on the mass-radius relationship and transmission and emission spectra. An interior-atmosphere model, combined with static structure calculations, tracks the evolving radius of a rocky mantle that is outgassing CO$_2$ and H$_2$O. Synthetic emission and transmission spectra are generated for CO$_2$ and H$_2$O dominated atmospheres. Atmospheres dominated by CO$_2$ suppress the outgassing of H$_2$O to a greater extent than previously realised, as previous studies have applied an erroneous relationship between volatile mass and partial pressure. We therefore predict more H$_2$O can be retained by the interior during the later stages of magma ocean crystallisation. Furthermore, formation of a lid at the surface can tie outgassing of H$_2$O to the efficiency of heat transport through the lid, rather than the atmosphere's radiative timescale. Contraction of the mantle as it solidifies gives $\sim5\%$ radius decrease, which can partly be offset by addition of a relatively light species to the atmosphere. We conclude that a molten silicate mantle can increase the radius of a terrestrial planet by around $5\%$ compared to its solid counterpart, or equivalently account for a $13\%$ decrease in bulk density. An outgassing atmosphere can perturb the total radius according to its speciation. Atmospheres of terrestrial planets around M-stars that are dominated by CO$_2$ or H$_2$O can be distinguished by observing facilities with extended wavelength coverage (e.g., JWST).

The energy balance of a partially molten rocky planet can be expressed as a non-linear diffusion equation using mixing length theory to quantify heat transport by both convection and mixing of the melt and solid phases. In this formulation the effective or eddy diffusivity depends on the entropy gradient, $\partial S/\partial r$, as well as entropy. First we present a simplified model with semi-analytical solutions, highlighting the large dynamic range of $\partial S/\partial r$, around 12 orders of magnitude, for physically-relevant parameters. It also elucidates the thermal structure of a magma ocean during the earliest stage of crystal formation. This motivates the development of a simple, stable numerical scheme able to capture the large dynamic range of $\partial S/\partial r$ and provide a flexible and robust method for time-integrating the energy equation. We then consider a full model including energy fluxes associated with convection, mixing, gravitational separation, and conduction that all depend on the thermophysical properties of the melt and solid phases. This model is discretised and evolved by applying the finite volume method (FVM), allowing for extended precision calculations and using $\partial S/\partial r$ as the solution variable. The FVM is well-suited to this problem since it is naturally energy conserving, flexible, and intuitive to incorporate arbitrary non-linear fluxes that rely on lookup data. Special attention is given to the numerically challenging scenario in which crystals first form in the centre of a magma ocean. Our computational framework is immediately applicable to modelling high melt fraction phenomena in Earth and planetary science research. Furthermore, it provides a template for solving similar non-linear diffusion equations arising in other disciplines, particularly for non-linear functional forms of the diffusion coefficient.