Compressible Convection Research Articles

Mass–Radius Relationships and Contraction of Condensed Planets by Cooling or Despinning

Condensed planets contract or expand as their temperature changes. With the exception of the effect of phase changes, this phenomenon is generally interpreted as being solely related to the thermal expansivity of the planet’s components. However, changes in density affect pressure and gravity and, consequently, the planet’s compressibility. A planet’s radius is also linked to its rate of rotation. Here again, changes in pressure, gravity, and compressibility are coupled. In this article we clarify how the radius of a condensed planet changes with temperature and rotation, using a simple and rigorous thermodynamic model. We consider condensed materials to obey a simple equation of state which generalizes a polytopic EoS as temperature varies. Using this equation, we build simple models of condensed planet’s interiors including exoplanets, derive their mass–radius relationships, and study the dependence of their radius on temperature and rotation rate. We show that it depends crucially on the value of ρ s gR/K s (ρ s being surface density, g gravity, R radius, K s surface incompressibility). This nondimensional number is also the ratio of the dissipation number which appears in compressible convection and the Gruneïsen mineralogic parameter. While the radius of small planets depends on temperature, this is not the case for large planets with large dissipation numbers; Earth and a super-Earth like CoRoT-7b are in something of an intermediate state, with a moderately temperature-dependent radius. Similarly, while the radius of these two planets is a function of their rotation rates, this is not the case for smaller or larger planets.

Internally heated and fully compressible convection: Flow morphology and scaling laws

In stars and planets natural processes heat convective flows in the bulk of a convective region rather than at hard boundaries. By characterizing how convective dynamics are determined by the strength of an internal heating source we can gain insight into the processes driving astrophysical convection. Internally heated convection has been studied extensively in incompressible fluids, but the effects of stratification and compressibility have not been examined in detail. In this work, we study fully compressible convection driven by a spatially uniform heating source in two-dimensional (2D) and 3D Cartesian, hydrodynamic simulations. We use a fixed temperature upper boundary condition which results in a system that is internally heated in the bulk and cooled at the top. We find that the flow speed, as measured by the Mach number, and turbulence, as measured by the Reynolds number, can be independently controlled by separately varying the characteristic temperature gradient from internal heating and the diffusivities. Two-dimensional simulations at a fixed Mach number (flow speed) demonstrate consistent power at low wave number as diffusivities are decreased. We observe convection where the velocity distribution is skewed towards cold, fast downflows and that the flow speed is related to the length scale and entropy gradient of the upper boundary where the downflows are driven. We additionally find a heat transport scaling law which is consistent with prior incompressible work. Published by the American Physical Society 2024

Convective scale and subadiabatic layers in simulations of rotating compressible convection

Context. Rotation is thought to influence the size of convective eddies and the efficiency of convective energy transport in the deep convection zones of stars. Rotationally constrained convection has been invoked to explain the lack of large-scale power in observations of solar flows. Aims. Our main aims are to quantify the effects of rotation on the scale of convective eddies and velocity as well as the depths of convective overshoot and subadiabatic Deardorff layers. Methods. We ran moderately turbulent three-dimensional hydrodynamic simulations of rotating convection in local Cartesian domains. The rotation rate and luminosity of the simulations were varied in order to probe the dependency of the results on Coriolis, Mach, and Richardson numbers measuring the influences of rotation, compressibility, and stiffness of the radiative layer. The results were compared with theoretical scaling results that assume a balance between Coriolis, inertial, and buoyancy (Archimedean) forces, also referred to as the CIA balance. Results. The horizontal scale of convective eddies decreases as rotation increases, and it ultimately reaches a rotationally constrained regime consistent with the CIA balance. Using a new measure of the rotational influence on the system, we found that even the deep parts of the solar convection zone are not in the rotationally constrained regime. The simulations captured the slowly and rapidly rotating scaling laws predicted by theory, and the Sun appears to be in between these two regimes. Both the overshooting depth and the extent of the Deardorff layer decrease as rotation becomes more rapid. For sufficiently rapid rotation, the Deardorff layer is absent due to the symmetrisation of upflows and downflows. However, for the most rapidly rotating cases, the overshooting increases again due to unrealistically large Richardson numbers that allow convective columns to penetrate deep into the radiative layer. Conclusions. Relating the simulations with the Sun suggests that the convective scale, even in the deep parts of the Sun, is only mildly affected by rotation and that some other mechanism is needed to explain the lack of strong large-scale flows in the Sun. Taking the current results at face value, the overshoot and Deardorff layers are estimated to span roughly 5% of the pressure scale height at the base of the convection zone in the Sun.

Strongly superadiabatic and stratified limits of compressible convection

Strongly superadiabatic and stratified limits of compressible convection

Compressible turbulent convection in highly stratified adiabatic background

Buoyancy-driven turbulent convection leads to a fully compressible flow with a prominent top-down asymmetry of first- and second-order statistics when the adiabatic equilibrium profiles of temperature, density and pressure change very strongly across the convection layer. The growth of this asymmetry and the formation of an increasingly thicker stabilized sublayer with a slightly negative mean convective heat flux$J_c(z)$at the top of the convection zone is reported here by a series of highly resolved three-dimensional direct numerical simulations beyond the Oberbeck–Boussinesq and anelastic limits for dimensionless dissipation numbers,$0.1 \le D\le 0.8$, at fixed Rayleigh number$Ra=10^6$and superadiabaticity$\epsilon =0.1$. The highly stratified compressible convection regime appears for$D > D_{crit}\approx 0.65$, when density fluctuations collapse to those of pressure; it is characterized by an up to nearly 50 % reduced global turbulent heat transfer and a sparse network of focused thin and sheet-like thermal plumes falling through the top sublayer deep into the bulk.

Compressible convection in super-earths

The radial density of planets increases with depth due to compressibility, leading to impacts on their convective dynamics. To account for these effects, including the presence of a quasi-adiabatic temperature profile and entropy sources due to dissipation, the compressibility is expressed through a dissipation number, D, proportional to the planet's radius and gravity. In Earth's mantle, compressibility effects are moderate, but in large rocky or liquid exoplanets (super-earths), the dissipation number can become very large. This paper explores the properties of compressible convection when the dissipation number is significant. We start by selecting a simple Murnaghan equation of state that embodies the fundamental properties of condensed matter at planetary conditions. Next, we analyze the characteristics of adiabatic profiles and demonstrate that the ratio between the bottom and top adiabatic temperatures is relatively small and probably less than 2. We examine the marginal stability of compressible mantles and reveal that they can undergo convection with either positive or negative superadiabatic Rayleigh numbers. Lastly, we delve into simulations of convection in 2D Cartesian geometry performed using the exact equations of mechanics, neglecting inertia (infinite Prandtl number case), and examine their consequences for super-earth dynamics.

A playground for compressible natural convection with a nearly uniform density

In the quest to understand the basic universal features of compressible convection, one would like to disentangle genuine consequences of compression from spatial variations of transport properties. For instance, one may choose to consider a fluid with uniform dynamic viscosity, but, then, compressible effects will generate a density gradient and consequently the kinematic viscosity will not be uniform. In the present work, we consider a very peculiar equation of state, whereby entropy is solely dependent on density, so that a nearly isentropic fluid domain is nearly isochoric. Within this class of equations of state, there is a thermal adiabatic gradient and a key property of compressible convection is still present, namely its capacity to viscously dissipate a large fraction of the thermal energy involved, of the order of the well-named dissipation number. In the anelastic approximation, under the assumption of an infinite Prandtl number, the number of governing parameters can be brought down to two, the Rayleigh number and the dissipation number. This framework is proposed as a playground for compressible convection, an opportunity to extend the vast corpus of theoretical analyses on the Oberbeck–Boussinesq equations regarding stability, bifurcations or the determination of upper bounds for the turbulent heat transfer. Here, in a two-dimensional geometry, we concentrate on the structure of numerical solutions. For all Rayleigh numbers, a change in the vertical temperature profile is observed in the range of dissipation number between $0$ and less than $0.4$ , associated with the weakening of ascending plumes. For larger dissipation numbers, the heat flux dependence on this number is found to be well predicted by Malkus's model of critical layers. For dissipation numbers of order unity, and large Rayleigh numbers, dissipation becomes related to the entropy heat flux at each depth, so that the vertical dissipation profile can be predicted, and so does the total ratio of dissipation to convective heat flux.

Fully compressible convection for planetary mantles

SUMMARY The numerical simulations of convection inside the mantle of the Earth or of terrestrial planets have been based on approximate equations of fluid dynamics. A common approximation is the neglect of the inertia term which is certainly reasonable as the Reynolds number of silicate mantles, or their inverse Prandtl number, are infinitesimally small. However various other simplifications are made which we discuss in this paper. The crudest approximation that can be done is the Boussinesq approximation (BA) where the various parameters are constant and the variations of density are only included in the buoyancy term and assumed to be proportional to temperature with a constant thermal expansivity. The variations of density with pressure and the related physical consequences (mostly the presence of an adiabatic temperature gradient and of dissipation) are usually accounted for by using an anelastic approximation (AA) initially developed for astrophysical and atmospheric situations. The BA and AA cases provide simplified but self-consistent systems of differential equations. Intermediate approximations are also common in the geophysical literature although they are invariably associated with theoretical inconsistencies (non-conservation of total energy, non-conservation of statistically steady state heat flow with depth, momentum and entropy equations implying inconsistent dissipations). We show that, in the infinite Prandtl number case, solving the fully compressible (FC) equations of convection with a realistic equation of state (EoS) is however not much more difficult or numerically challenging than solving the approximate cases. We compare various statistical properties of the Boussinesq, AA and FC simulations in 2-D simulations. We point to an inconsistency of the AA approximation when the two heat capacities are assumed constant. We suggest that at high Rayleigh number, the profile of dissipation in a convective mantle can be directly related to the surface heat flux. Our results are mostly discussed in the framework of mantle convection but the EoS we used is flexible enough to be applied for convection in icy planets or in the inner core.

An Interpolated Bounce Back Thermable Method for Simulating Solid Particles Dynamics in a Viscous Medium

In this paper we discuss the mathematical and computer modeling of non-isothermal two-phase flows with suspended particles. Natural convection between an outer cubical cavity and an inner hot sphere is investigated. To simulate heat fluxes loaded with particles, a thermal model of the lattice Boltzmann equation in combination with the interpolated bounce back method (TLBM-IBB) has been developed. In TLBM-IBB, IBB is used to process liquid-solid interfaces, and TLBM is used to simulate the heat flow of a fluid. The momentum exchange method is used to calculate the hydrodynamic force on the particle surface. Simulation performed for a range of Rayleigh numbers . The accuracy and efficiency of the existing method is demonstrated by the example of solving the test problem of natural convection around a stationary particle and three-dimensional compressible natural convection in a square cavity filled with air, which has a hot wall on the left and a cold wall on the right, and two horizontal walls are adiabatic. The results obtained are in good agreement with the experimental and numerical results of other authors.

Prandtl number dependence of stellar convection: Flow statistics and convective energy transport

Context. The ratio of kinematic viscosity to thermal diffusivity, the Prandtl number, is much smaller than unity in stellar convection zones. Aims. The main goal of this work is to study the statistics of convective flows and energy transport as functions of the Prandtl number. Methods. Three-dimensional numerical simulations of compressible non-rotating hydrodynamic convection in Cartesian geometry are used. The convection zone (CZ) is embedded between two stably stratified layers. The dominant contribution to the diffusion of entropy fluctuations comes in most cases from a subgrid-scale diffusivity whereas the mean radiative energy flux is mediated by a diffusive flux employing Kramers opacity law. Here, we study the statistics and transport properties of up- and downflows separately. Results. The volume-averaged rms velocity increases with decreasing Prandtl number. At the same time, the filling factor of downflows decreases and leads to, on average, stronger downflows at lower Prandtl numbers. This results in a strong dependence of convective overshooting on the Prandtl number. Velocity power spectra do not show marked changes as a function of Prandtl number except near the base of the convective layer where the dominance of vertical flows is more pronounced. At the highest Reynolds numbers, the velocity power spectra are more compatible with the Bolgiano-Obukhov k−11/5 than the Kolmogorov-Obukhov k−5/3 scaling. The horizontally averaged convected energy flux (F̅conv), which is the sum of the enthalpy (F̅enth) and kinetic energy fluxes (F̅kin), is independent of the Prandtl number within the CZ. However, the absolute values of F̅enth and F̅kin increase monotonically with decreasing Prandtl number. Furthermore, F̅enth and F̅kin have opposite signs for downflows and their sum F̅↓conv diminishes with Prandtl number. Thus, the upflows (downflows) are the dominant contribution to the convected flux at low (high) Prandtl numbers. These results are similar to those from Rayleigh-Benárd convection in the low Prandtl number regime where convection is vigorously turbulent but inefficient at transporting energy. Conclusions. The current results indicate a strong dependence of convective overshooting and energy flux on the Prandtl number. Numerical simulations of astrophysical convection often use a Prandtl number of unity because it is numerically convenient. The current results suggest that this can lead to misleading results and that the astrophysically relevant low Prandtl number regime is qualitatively different from the parameter regimes explored in typical contemporary simulations.

Radius ratio dependency of the instability of fully compressible convection in rapidly rotating spherical shells

Based on the fully compressible Navier–Stokes equations, the linear stability of thermal convection in rapidly rotating spherical shells of various radius ratios $\eta$ is studied for a wide range of Taylor number $Ta$, Prandtl number $Pr$ and the number of density scale height $N_\rho$. Besides the classical inertial mode and columnar mode, which are widely studied by the Boussinesq approximation and anelastic approximation, the quasi-geostrophic compressible mode is also identified in a wide range of $N_\rho$ and $Pr$ for all $\eta$ considered, and this mode mainly occurs in the convection with relatively small $Pr$ and large $N_\rho$. The instability processes are classified into five categories. In general, for the specified wavenumber $m$, the parameter space ($Pr, N_\rho$) of the fifth category, in which the base state loses stability via the quasi-geostrophic compressible mode and remains unstable, shrinks as $\eta$ increases. The asymptotic scaling behaviours of the critical Rayleigh numbers $Ra_c$ and corresponding wavenumbers $m_c$ to $Ta$ are found at different $\eta$ for the same instability mode. As $\eta$ increases, the flow stability is strengthened. Furthermore, the linearized perturbation equations and Reynolds–Orr equation are employed to quantitatively analyse the mechanical mechanisms and flow instability mechanisms of different modes. In the quasi-geostrophic compressible mode, the time-derivative term of disturbance density in the continuity equation and the diffusion term of disturbance temperature in the energy equation are found to be critical, while in the columnar and inertial modes, they can generally be ignored. Because the time-derivative term of the disturbance density in the continuity equation cannot be ignored, the anelastic approximation fails to capture the instability mode in the small-$Pr$ and large-$N_\rho$ system, where convection onset is dominated by the quasi-geostrophic compressible mode. However, all the modes are primarily governed by the balance between the Coriolis force and the pressure gradient, based on the momentum equation. Physically, the most important difference between the quasi-geostrophic compressible mode and the columnar mode is the role played by the disturbance pressure. The disturbance pressure performs negative work for the former mode, which appears to stabilize the flow, while it destabilizes the flow for the latter mode. As $\eta$ increases, in the former mode the relative work performed by the disturbance pressure increases and in the latter mode decreases.

Development of unsteady natural convection in a square cavity under large temperature difference

To investigate how the nonuniform fluid density distribution caused by large temperature variations affects the development of unsteady natural convection, we perform a series of direct numerical simulations of two-dimensional compressible natural convection in an air-filled square cavity. The cavity has a hot wall on the left and a cold wall on the right, and two horizontal walls are adiabatic. The simulations are done using a kinetic approach based on a modeled Boltzmann equation, from which the fully compressible Navier–Stokes–Fourier equations are recovered. No Boussinesq approximation or low Mach number approximation is made. An extra source term is introduced to adjust the fluid Prandtl number. Simulations are performed for a range of Rayleigh numbers (107−109) with a fixed dimensionless temperature difference of ε=0.6 to determine the critical Rayleigh number and study the development of unsteady flow. To illustrate the instability mechanism, instantaneous fluctuation field, time trace of temperature, and velocity at selected monitoring points, the spectrum and other statistics are presented and discussed. As expected, significant differences are observed between the instability of compressible natural convection and the Boussinesq-type natural convection. With a large temperature difference, the transition to unsteady flow is asymmetric for the flows near the hot wall and cold wall. For the Rayleigh number range we studied, the cold wall region is dominated by low-frequency impact instability of the boundary thermal jet at the bottom corner. For the hot wall region, besides the upper corner impact instability, a boundary layer instability featuring high-frequency oscillations is observed.

Designing a consistent implementation of the discrete unified gas-kinetic scheme for the simulation of three-dimensional compressible natural convection

Discrete unified gas-kinetic scheme (DUGKS) has been developed as a robust and accurate approach for thermal compressible flow simulations; however, designing an efficient and accurate lattice velocity model to take full advantage of DUGKS remains a challenge. In this study, we apply DUGKS to simulate three-dimensional compressible natural convection in an enclosure with a large temperature difference, without making the Boussinesq approximation. The Chapman–Enskog analysis indicates that the fourth-order moments of equilibrium is needed for the heat flux evaluation in the energy equation, implying that the fourth-order Hermite expansion of equilibrium and thus at least an eighth-order Gauss–Hermite quadrature are needed for accurate simulation of the Navier–Stokes–Fourier system. For this purpose, a highly efficient lattice velocity model, D3Q77A9, is derived, which provides a Gauss–Hermite quadrature of ninth-order accuracy in three dimensions. The accuracy of this D3Q77A9 model is demonstrated by simulating compressible natural convection flows in both two-dimensional and three-dimensional cavities. An error analysis is performed to emphasize the importance of combining a quadrature with an adequate degree of precision and a proper order of Hermite expansion of the equilibrium distribution.

An improved discrete unified gas kinetic scheme for simulating compressible natural convection flows

Discrete unified gas-kinetic scheme (DUGKS) has been developed recently as a general method for simulating flows at all Knudsen numbers. In this study, we extend DUGKS to simulate fully compressible thermal flows. We introduce a source term to the Boltzmann equation with the Bhatnagar-Gross-Krook (BGK) collision model [1] to adjust heat flux and thus the Prandtl number. The fully compressible Navier-Stokes equations can be recovered by the current model. As a mesoscopic CFD approach, it requires an accurate mesoscopic implementation of the boundary conditions. Using the Chapman-Enskog approximation, we derive the “bounce-back” expressions for both temperature and velocity distribution functions, which reveal the need to consider coupling terms between the velocity and thermal fields. To validate our scheme, we first reproduce the Boussinesq flow results by simulating natural convection in a square cavity with a small temperature difference (ϵ=0.01) and a low Mach number. Then we perform simulations of steady natural convection (Ra=1.0×106) in a square cavity with differentially heated side walls and a large temperature difference (ϵ=0.6), where the Boussinesq approximation becomes invalid. Temperature, velocity profiles, and Nusselt number distribution are obtained and compared with the benchmark results from the literature. Finally, the unsteady compressible natural convection with Ra=5.0×109,ϵ=0.6 is studied and the turbulent fluctuation statistics are computed and analyzed.

The effect of temperature-dependent viscosity and thermal conductivity on the onset of compressible convection

The linear equations of thermal convection in a compressible fluid with non-constant transport coefficients are derived. The criterion for the onset of convection is established, based on linear stability analysis, for a range of different temperature-dependent profiles of thermal conductivity and viscosity. Temperature-dependent transport coefficients are shown to lead to a more complex behaviour than their constant counterparts, and modifies the stability condition of the fluid. When the Rayleigh number is defined in terms of the mid-layer physical properties and the temperature gradient at the top is held constant, increasing the temperature-dependence of thermal conductivity is found to raise the critical Rayleigh number dramatically, as the convective disturbance is then concentrated mainly at the top of the layer. In contrast, for viscosity a more subtle effect on stability is identified.

Comment on “Compressible Coupled Natural Convection and Volumetric Radiation in a Square Cavity”

Comment on “Compressible Coupled Natural Convection and Volumetric Radiation in a Square Cavity”

Comparison of 2D and 3D compressible convection in a pre-main sequence star

Context.A 1D description of stellar dynamics is at the basis of stellar evolution modeling. Designed to investigate open problems in stellar evolution, the MUltidimensional Stellar Implicit Code expands a realistic 1D profile of a star’s internal structure to examine the interior dynamics of a specific star through either 2D or 3D hydrodynamic simulations.Aims.Extending our recent studies of 2D stellar convection to 3D stellar convection, we aim to compare 3D hydrodynamic simulations to identically set-up 2D simulations, for a realistic pre-main sequence star.Methods.We compare statistical quantities related to convective flows including: average velocity, vorticity, local enstrophy, and penetration depth beneath a convection zone. These statistics were produced during stationary, steady-state compressible convection in the star’s convection zone.Results.Our simulations confirm the common result that 2D simulations of stellar convection have a higher magnitude of velocity on average than 3D simulations. Boundary conditions and the extent of the spherical shell can affect the magnitude and variability of convective velocities. The difference between 2D and 3D velocities is dependent on these background points; in our simulations this can have an effect as large as the difference resulting from the dimensionality of the simulation. Nevertheless, radial velocities near the convective boundary are comparable in our 2D and 3D simulations. The average local enstrophy of the flow is lower for 2D simulations than for 3D simulations, indicating a different shape and structuring of 3D stellar convection. We performed a statistical analysis of the depth of convective penetration below the convection zone using the model proposed in our recent study (Pratt et al. 2017, A&A, 604, A125). That statistical model was developed based on 2D simulations, which allowed us to examine longer times and higher radial resolution than are possible in 3D. Here, we analyze the convective penetration in 3D simulations, and compare the results to identically set-up 2D simulations. In 3D simulations, the penetration depth is as large as the penetration depth calculated from 2D simulations.

Upward Overshooting in Turbulent Compressible Convection. II. Simulations at Large Relative Stability Parameters

Abstract In this paper, we present the results of three-dimensional numerical simulation of upward overshooting in turbulent compressible convection at large relative stability parameter S. Similar to the previous simulations at small S, we find that the convectively stable zone can be partitioned into three layers: the thermal adjustment layer, the turbulent dissipation layer, and the thermal dissipation layer. Despite of this similarity, there exist significant differences in several aspects. First, for small S, the thermal structure is altered considerably near the interface between the convectively unstable and stable zones. For extremely large S, the thermal structure is only slightly changed. Second, the overshooting distance decreases at small S, but it can increase when S is large enough. Third, for small S, the fluid motions tend to be less active when S increases. However, the fluid motions can be more active when S is large enough. We find that the structure of counter cells has a significant impact on the penetration depth.

Upward Overshooting in Turbulent Compressible Convection. III. Calibrate Parameters for One-dimensional Reynolds Stress Model

Abstract In this paper, we calibrate the coefficients for the one-dimensional Reynolds stress model (RSM) with the data generated from the three-dimensional (3D) numerical simulations of upward overshooting in turbulent compressible convection. It has been found that the calibrated convective and isotropic coefficients are almost the same as those calibrated in the pure convection zone. However, the calibrated diffusive coefficients differ significantly from those calibrated in the pure convection zone. We suspect that the diffusive effect induced by the boundary is stronger than by the adjacent stable zone. We have checked the validity of the downgradient approximation (DGA). We find that the prediction of the DGA on the third-order moments (TOMs) is unsatisfactory. However, the prediction on their derivatives is much better. It explains why the performance of the RSM is reasonable in application to the real stars. With the calibrated coefficients, we have solved the full set of nonlocal turbulent equations on RSM. We find that the RSM has successfully produced the thermal adjustment layer and turbulent dissipation layer, which were identified in the 3D numerical simulations. We suggest to use the inflection point of the auto-correlation of temperature perturbation and the Péclet number as the indicators on measuring the extents of the thermal adjustment layer and turbulent dissipation layer, respectively. This result may offer a practical guidance on the application of the RSM in 1D stellar structure and evolution models.

On formulations of compressible mantle convection

SUMMARY Mantle convection and long-term lithosphere dynamics in the Earth and other planets can be treated as the slow deformation of a highly viscous fluid, and as such can be described using the compressible Navier–Stokes equations. Since on Earth-sized planets the influence of compressibility is not a dominant effect, density deviations from a reference profile are at most on the order of a few percent and using the full governing equations poses numerical challenges, most modelling studies have simplified the governing equations. Common approximations assume a temporally constant, but depth-dependent reference profile for the density (the anelastic liquid approximation), or drop compressibility altogether and use a constant reference density (the Boussinesq approximation). In most previous studies of mantle convection and crustal dynamics, one can assume that the error introduced by these approximations was small compared to the errors that resulted from poorly constrained material behaviour and limited numerical accuracy. However, as model parametrizations have become more realistic, and model resolution has improved, this may no longer be the case and the error due to using simplified conservation equations might no longer be negligible: while such approximations may be reasonable for models of mantle plumes or slabs traversing the whole mantle, they may be unsatisfactory for layered materials experiencing phase transitions or materials undergoing significant heating or cooling. For example, at boundary layers or close to dynamically changing density gradients, the error arising from the use of the aforementioned compressibility approximations can be the dominant error source, and common approximations may fail to capture the physical behaviour of interest. In this paper, we discuss new formulations of the continuity equation that include dynamic density variations due to temperature, pressure and composition without using a reference profile for the density. We quantify the improvement in accuracy relative to existing formulations in a number of benchmark models and evaluate for which practical applications these effects are important. Finally, we consider numerical aspects of the new formulations. We implement and test these formulations in the freely available community software aspect, and use this code for our numerical experiments.