15
$\begingroup$

I'm thinking of a planet in a Goldilocks zone similar to Earth's, with a similar atmosphere, and similar atmospheric pressure and temperature at the surface. Gravity would be variable, based on the mass needed to sustain the kinds of pressures to form exotic ices at the core.

I realize there are at least a couple of similar questions -

Could a planet made completely of water exist?

What would happen at the core of a water world?

  • but I'm specifically wondering about the necessary mass to achieve these matter states at the core.

Thanks!

$\endgroup$
6
  • $\begingroup$ worldbuilding.stackexchange.com/a/54139/8068 Maybe solid-matter physicists would know (or at least have the formulas). A direct questionis, "How much water is required to create gravitational pressure of 5GPa (aka 725000 psi)?" $\endgroup$
    – RonJohn
    Mar 11, 2018 at 6:16
  • $\begingroup$ Thanks. I'll consider cross-posting it, depending on the reception here. Maybe sometime in the next week. $\endgroup$ Mar 11, 2018 at 7:20
  • $\begingroup$ @RonJohn, I'm not sure where you're getting the 5 GPa from. From the phase diagram, that looks like the ice VIII -> ice VII transition at 0C, but you're not going to have straight water on top of your ice VII at that temperature -- you're going to have layers of ice VI and ice V (and a layer of ice 1h on top of the water, but that one is thin enough to ignore). $\endgroup$
    – Mark
    Mar 11, 2018 at 8:01
  • 1
    $\begingroup$ @Mark Looking at the phase diagram, it will also depend on the temperature at the core too. So knowing both the temperature & pressure of the core will be critical or at least calculating the column of water necessary to create those conditions. From that, estimate the mass of the water planet. $\endgroup$
    – a4android
    Mar 11, 2018 at 12:05
  • $\begingroup$ @Mark by following one of the links in the article. (I knew I should have quoted it...) I'm not in the mood to track it down, though. $\endgroup$
    – RonJohn
    Mar 11, 2018 at 15:32

3 Answers 3

12
$\begingroup$

Your question is similar to another one I answered, and I'll be borrowing heavily from that.

Ice VII

I can start off the answers with some simplifying assumptions, but someone else may have a better idea of the intricacies involved in this calculation. The specific assumptions I'll make are:

  1. Constant temperature

  2. Water is incompressible, and on planetary scales has an average density of ~1.5g/cm$^3$ (see my other answer for justification)

With these assumptions, this really just becomes a plugging-in-the-numbers problem.

Here's the water phase diagram I'll use to talk about the rest of this problem:

Water phase diagram

Given our first assumption, I'm going to choose a balmy 350K for the planet's water and ice. In the diagram above, we can see that ice VII shows up at ~2GPa. So the question becomes, how much water does it take to get a 2GPa pressure at the core?

Well, given our second assumption, it would take 200 kilometers of water to reach 2GPa given the classic conversion 101kPa/10m. With that info, we can calculate the mass of a the planet via the equation

$$m = density*volume = \rho*\frac{4\pi r^3}{3}$$

$$(1.5\frac{g}{cm^3}*10^{15}\frac{cm^3}{km^3})*(\frac{4\pi 200^3\ km^3}{3}) = 5*10^{22}g$$

$5*10^{19}$ kg

Cool! From here, we can see that's in the range of large asteroids or small moons.

Of course, this is fundamentally a ballpark estimate, but I'd say it's accurate to about an order of magnitude. As pointed out in the comments on my other answer, the constant temperature assumption is possibly quite valid depending on how your planet formed and how old it is. The "average" water density leaves a lot of room for error, but I wasn't confident enough in my calculus to do the full derivative (change in bulk modulus with respect to pressure as depth increases).

Ice X

This one gets a lot more complicated because we're working with two shells here, but we have a similar set of assumptions:

  1. Constant temperature

  2. Water is incompressible, and on planetary scales has an average density of ~1.5g/cm$^3$

  3. Ice VII is incompressible, and on planetary scales has an average density of ~2.3g/cm$^3$

At 350K, we reach ice X at about 50 GPa and can answer the question similarly to the one above- how thick does the ice VII need to be to reach this pressure? We already know that we'll have 200km liquid water along the surface, so the core is the only new thing here.

To get an additional 48 GPa due to Ice VII, we need approximately 2000 additional kilometers:

$$h = \frac{48*10^9}{2300*9.8} = 2130km$$

Take this estimate with a large grain of salt- $g$ wouldn't be constant through a planet's core, but it would rather depend on the mass of the planet and the distance from the surface, which means we're into those nasty differential equations again. Man, no wonder physicists are angry all the time.

With this, we can again calculate our mass with the equation:

$$m = (V_{core}*\rho_{core}+V_{ocean}*\rho_{ocean})$$

which, if we plug in our numbers and solve it correctly, returns

$8.9*10^{22}$ kg

Which is about as large as the biggest moons and in the range of the smallest planets. Good question!

$\endgroup$
4
  • $\begingroup$ Thank you! This is actually for a fantasy setting, but I want as little as possible to be hand waved with magic. $\endgroup$ Mar 12, 2018 at 22:39
  • $\begingroup$ Either I messed something up in my answer or you can't ignore the variation of gravity and pressure with depth: I get masses that are two to three orders of magnitude higher than yours, and radiuses significantly higher. $\endgroup$
    – Mark
    Mar 13, 2018 at 5:56
  • 2
    $\begingroup$ I've done some work and I came up with results that agree with @Mark's fairly well, and I think he's right - the variation of gravity at different radii seems like the only possible culprit. $\endgroup$
    – HDE 226868
    May 9, 2019 at 14:05
  • $\begingroup$ I kinda feel like I should upvote only the two "correctish" answers. But I also feel compelled to upvote "good sportsmanship in answering". This does not help correct the relative ordering, but feels more "right". $\endgroup$ May 10, 2019 at 21:19
14
$\begingroup$

I decided to write a program to calculate this. It iteratively builds a planet from the core outwards in one-meter layers, calculating the gravity and adjusting the density of the current material for the pressure at each step.

#include <math.h>
#include <stdio.h>

/* Units are meters, kilograms, seconds */
const double G = 0.00000000006674;

struct
{
    double baseDensity;
    double bulkModulus;
    double lowerPressure;
    double upperPressure;
    const char *name;
} Properties[] = {
    {1000, 2200000000, 50000, 2000000000, "water"},
    {1500, 23900000000 , 2000000000, 50000000000, "ice VII"},
    {2500, 23900000000, 50000000000, 400000000000, "ice X"},
    {3000, 10000000000000000, 400000000000, 1000000000000, "ice XI"}    /* Mostly made-up, but it doesn't matter, because we've only got a 1-meter sphere of it. */
};

/* Calculate from the inside out. */
void CalculatePlanet(double *radius, double *mass)
{
    int currentMaterial = 3;    /* Start with a 1-meter layer of ice XI */
    double pressureNeeded = Properties[currentMaterial].lowerPressure;  /* We need to stack up material to produce this much pressure */
    *radius = 0;
    *mass = 0;

    while(currentMaterial >= 0)
    {
        double shellMass;   /* Mass of the shell */
        double shellPressure;   /* Pressure provided by the shell */
        double newRadius = *radius + 1.0;
        double shellDensity = (pressureNeeded * Properties[currentMaterial].baseDensity) / Properties[currentMaterial].bulkModulus + Properties[currentMaterial].baseDensity;

        /* Add a one-meter layer to the planet */
        shellMass = (((newRadius) * (newRadius) * (newRadius)) - (*radius * *radius * *radius)) * (4.0/3.0) * M_PI * shellDensity;

        shellPressure = G * *mass / (*radius * *radius) * shellDensity;
        if(isnan(shellPressure)) shellPressure = 0;

        pressureNeeded -= shellPressure;
        *mass += shellMass;
        *radius += 1.0;

        if(pressureNeeded < Properties[currentMaterial].lowerPressure)
        {
            printf("Layer: %i %lf %lf %lf %lf %lf\n", currentMaterial, shellPressure, pressureNeeded, shellMass, *mass, *radius);
            currentMaterial--;
        }
    }
}


int main(void)
{
    double mass = 0;
    double radius = 0;
    CalculatePlanet(&radius, &mass);

    double volume = radius * radius * radius * M_PI * 4.0 / 3.0;
    double density = mass / volume;
    double surfaceGravity = G * mass /(radius * radius);

    printf("Planet calculated.  Radius %.0lf meters, mass %.0lf kg, density %0lf kg/m3, gravity %lf m/s2\n", radius, mass, density, surfaceGravity);
}

Using the same 350K planet, presumed bulk moduli, and phase diagram as Dubukay, I get the following planets:

Core of water (done as a sanity check): radius 1 meters, mass 4189 kg, density $1000 kg/m^3$

Core of ice VII, surrounded by 2555498 meters of water: radius 2555499 meters, mass $8.98 * 10^{22}$ kg, density $1285 kg/m^3$, surface gravity $0.92 m/s^2$. Roughly the diameter of Mercury, but only a quarter as heavy.

Core of ice X, surrounded by 6013480 meters of ice VII and 349831 meters of water: radius 6363312 meters, mass $2.44 * 10^{24} kg$, density $2261 kg/m^3$, surface gravity $4.02 m/s^2$. About as large as Earth, but only 40% the mass.

Core of ice XI, surrounded by 2209965 meters of ice X, 2675055 meters of ice VII, and 301287 meters of water: radius 5186308 meters, mass $1.85 * 10^{24} kg$, density $3174 kg/m^3$, surface gravity $4.60 m/s^2$. A bit smaller than Earth, and only a third the mass.

Note that the planet with a core of ice X is larger than the planet with a core of ice XI. This isn't an error: ice X is far denser than ice VII; the reduced radius increases the gravity at all levels, making for higher pressures and densities.

$\endgroup$
2
  • $\begingroup$ Nice, +1! Shouldn't your last remark meaning that the planet with Ice X core would be unstable --- prone to collapse to the Ice XI-cored one? What a wave... ;-) (but well, I see that temperature ranges are vastly different). $\endgroup$
    – Rmano
    Nov 27, 2018 at 13:08
  • $\begingroup$ @Rmano, it's not unstable per se, but during planetary accretion, once you get enough material to form Ice X, you're going to get a fairly catastrophic transition. $\endgroup$
    – Mark
    Nov 28, 2018 at 4:42
9
+100
$\begingroup$

Summary

It turns out that even relatively low-mass ocean planets are capable of forming some of the exotic ices you name in their cores. Ice VII appears to form at the centers of planets of $0.015M_{\oplus}$ (Earth masses), while ice X forms at the centers of planets of $1.256M_{\oplus}$. Interestingly, despite the increase in mass by two orders of magnitude and the increase in central pressure by a factor of 25, these worlds have radii differing by only a factor of four. While there may be a temperature dependence, given the relative simplicity of water's phase diagram at $\sim300\text{ K}$, I suspect this should not be an issue, and the relevant equations of state are not temperature-dependent.

Theory

Since we have two competing answers (Dubukay's and Mark's) with vastly different results, I thought I would add a third method to see if I could come up with something similar. I went to Seager et al. 2008, my favorite set of models of the interiors of terrestrial exoplanets. Their setup assumes that the bodies are isothermal at low pressures - as Dubukay did - and uses equations of state of the form $$\rho(P)=\rho_0+cP^n\tag{11}$$ where $\rho$ is density, $P$ is pressure and $c$ and $n$ are composition-dependent constants; $n\approx0.5$ for most terrestrial worlds, but it does differ, which is important. This equation is essentially a modified polytrope, with the one major change being that $\rho(0)\neq0$, which would be true in a classic polytrope. For a pure $\text{H}_2\text{O}$ planet, $n=0.513$ and $c=0.00311$. When using these constants, bear in mind that pressure is in pascals, and density is in kilograms per cubic meter.

Seager et al. derive the following mass-radius relationship (I have numbered the equations as they are numbered in the paper): $$M(R)=\frac{4\pi}{3}R^3\left[\rho(P_c)-\frac{2}{5}\pi GR^2\rho_0^2f'(P_c)\right]\tag{31}$$ where $f(P)=cP^n$ and $P_c$ is the central pressure. It can be shown via hydrostatic equilibrium that $$P_c=\frac{3G}{8\pi}\frac{M^2}{R^4}\tag{27}$$ Given a desired central pressure, I can test various radii and corresponding masses and find the values I need.

We can check these results a different way: by numerical integration. The structure of any planet is governed by two key equations: $$\frac{dP}{dr}=-\frac{Gm\rho}{r^2}$$ $$\frac{dm}{dr}=4\pi r^2\rho$$ These are the equations of hydrostatic equilibrium and mass continuity. $r$ is a radial coordinate, measured from the center of the planet, and $m$ is the mass enclosed within $r$. By modeling the planet as a collection of progressively larger shells, and knowing the value of $P$ and $m$ in any given shell, we can find the value of $P$ and $m$ in the next shell via the Euler method: finding the change in these variables by multiplying their derivatives at a point by some step size $\Delta r$. This is essentially what Mark did, I think. I'm simply using a particular equation of state, rather than a bulk modulus.

Code

I wrote some fairly simple code in Python 3 to accomplish this. It only requires NumPy (as well as Matplotlib for auxiliary plots).

import numpy as np

earthMass = 5.97*10**(24) # kg
earthRadius = 6.371*10**(6) # m
G = 6.67*10**(-11) # gravitational constant, SI units

def rho(P,rho0,c,n):
    """Polytropic equation of state"""
    rho = rho0 + c*(P**n)
    return rho

def fprime(P,c,n):
    """Derivative of the first order contribution
    to the polytropic equation of state"""
    fprime = c*n*(P**(n-1))
    return fprime

def mass(R,rho0,c,n):
    """Compute planetary mass for a particular radius,
    given equation of state parameters for a particular
    composition."""
    Rscaled = R*earthRadius # convert to SI units
    Pc = (2*np.pi/3)*G*(Rscaled**2)*(rho0**2) # central pressure
    rho_mean = rho(Pc,rho0,c,n) - (2*np.pi/5)*G*(Rscaled**2)*(rho0**2)*fprime(Pc,c,n) # mean density
    Mscaled = (4*np.pi/3)*(Rscaled**3)*rho_mean
    Mp = Mscaled/earthMass # convert to Earth masses
    return Mp

def pressure(R,rho0,c,n):
    """Compute central pressure if radius is known"""
    M = mass(R,rho0,c,n)
    M = M*earthMass # convert to SI units
    R = R*earthRadius # convert to SI units
    Pc = (3*G/(8*np.pi))*(M**2)/(R**4)
    return Pc

def minimumMass(P,rho0,c,n):
    """Compute mass at which a particular central
    pressure is reached"""
    radii = np.logspace(-1,1,1000) # reasonable radius range
    i = 0
    r = radii[i]
    while pressure(r,rho0,c,n) < P:
        # Brute force check of various radii
        i += 1
        r = radii[i]
    return(mass(r,rho0,c,n))

def radius(M,rho0,c,n):
    """Compute radius which yields a given mass"""
    radii = np.logspace(-1,1,1000)
    i = 0
    r = radii[i]
    while mass(r,rho0,c,n) < M:
        # Brute force check of various radii
        i += 1
        r = radii[i]
    return r

pressureList = [2,50] # central pressures to check, in GPa

for p in pressureList:
    print('Central pressure: '+str(p)+' GPa.')
    print('  The required mass is '\
          +str('%.3f'%minimumMass(p*10**9,1460,0.00311,0.513))+\
          ' Earth masses.')
    print('  The required radius is '+\
          str('%.3f'%radius(minimumMass(p*10**9,1460,0.00311,\
              0.513),1460,0.00311,0.513))+' Earth radii.')

Here is my numerical integration code. It's written specifically for water worlds, so the equation of state parameters are not function arguments. If you want to, it can be generalized easily enough for any composition.

import numpy as np

earthMass = 5.97*10**(24) # kg
earthRadius = 6.371*10**(6) # m
G = 6.67*10**(-11) # gravitational constant, SI units

rho0 = 1460
c = 0.00311
n = 0.513

def dP(M,R,P,dR):
    """Compute change in pressure via hydrostatic
    equilibrium"""
    rho = rho0 + c*(P**n) # density
    dP = -((G*M*rho)/(R**2))*dR
    return dP

def dM(R,P,dR):
    """Compute change in mass via mass continuity
    equation"""
    rho = rho0 + c*(P**n) # density
    dM = 4*np.pi*(R**2)*rho*dR
    return dM

def integrator(Pc,dR):
    """Numerically integrate differential equations
    to construct the planet"""
    P = [Pc,Pc]
    M = [0,0]
    R = [0,dR]
    # To avoid singularities at r = 0, I really
    # start the code at one step, r = dR. I assume
    # that this step is small enough that the mass
    # and pressure don't change significantly.

    while P[-1] > 0:
        # The surface of the planet is where P = 0
        m = M[-1]
        r = R[-1]
        p = P[-1]
        deltaR = 1
        deltaP = dP(m,r,p,deltaR)
        deltaM = dM(r,p,deltaR)
        P.append(P[-1]+deltaP)
        M.append(M[-1]+deltaM)
        R.append(R[-1]+deltaR)

    return M, R, P

pressureList = [2,50] # central pressures to check, in GPa

for p in pressureList:
    massList, radiusList, pressureList = integrator(p*(10**9),1)
    M = massList[-1]/earthMass
    R = radiusList[-1]/earthRadius
    print('Central pressure: '+str(p)+' GPa.')
    print('  The required mass is '+str('%.3f'%M)+\
          ' Earth masses.')
    print('  The required radius is '+str('%.3f'%R)+\
          ' Earth radii.')

Results

I chose a central pressure of $P_c=2\text{ GPa}$ for ice VII and $P_c=50\text{ GPa}$ for ice X, as Dubukay and Mark did. For both cases, my results agreed with Mark's to within an order of magnitude; the discrepancy with Dubukay's numbers still remains:

$$\begin{array}{|c|c|c|c|} \hline \text{} & \text{Ice VII} & \text{Ice X}\\ \hline \text{Dubukay} & M=8.327\times10^{-6}M_{\oplus} & M=0.0149M_{\oplus} \newline & R=0.0313R_{\oplus} & R=0.334R_{\oplus}\\ \hline \text{Mark} & M=0.0149M_{\oplus} & M=0.409M_{\oplus} \newline & R=0.401R_{\oplus} & R=0.998R_{\oplus}\\ \hline \text{Analytical} & M=0.0154M_{\oplus} & M=1.256M_{\oplus} \newline \text{models} & R=0.377R_{\oplus} & R=1.525R_{\oplus}\\ \hline \text{Numerical} & M=0.015M_{\oplus} & M=0.959M_{\oplus} \newline \text{integration} & R=0.372R_{\oplus} & R=1.389R_{\oplus}\\ \hline \end{array}$$

Both of my ice VII models agree very closely with Mark's, and my ice X models are only off by a factor of a few. The numerical integration does not match the analytical models, which worries me a little bit, but the discrepancy is not overly serious, and I'll do some poking around to see if I can find the problem. I'm happy enough to get within an order of magnitude in astronomy, so I'll consider all of this a victory. Here's a plot of my analytical results, with the terrestrial planets of the Solar System for comparison, as well as a curve of silicate planets ($\text{MgSiO}_3$):

Plot showing our two planets, as well as Mercury, Venus, Earth and Mars

What's going on?

This does shed some light on the different answers because a more detailed look at the theory rules out possible reasons for the discrepancy. The equations of state I used are isothermal; the other answers assume the same. Similarly, simple plots of density within these planets indicate that the weak dependence on pressure indeed justifies Dubukay's assumption of incompressibility. Both cases see perhaps a 10% change in density from the inner core to the surface - hardly enough to cause a discrepancy of three orders of magnitude. Indeed, at these pressures, most worlds should be quite incompressible.

I suspect that the key problem with Dubukay's answer is the assumption that the pressure-depth relationship doesn't change based on depth - and it likely does. By plotting the density inside each planet, we can see that it changes only slightly for the ice VII planet and a bit more for the ice X planet:

Plot of density profile of the planets

Now, the gravitational acceleration $g(r)$ at a radius $r$ scales as $g\propto\bar{\rho}r$, where $\bar{\rho}$ is the mean density inside $r$. The deviations from constant density are small for most regions onside the planet, so we should expect $g(r)$ to be fairly linear, and it is (closer to linear for the ice VII planet, which has a more uniform density profile):

Plot of gravitational acceleration inside the planets

Therefore, the simple depth-to-pressure conversion is inaccurate far from the surface. I also suspect that the core-ocean model is a little too simple.

$\endgroup$
1
  • 2
    $\begingroup$ Dude, nice. I’ll bounty you and Mark some of the rep I gained from this answer since my assumptions are too far off to be considered “correct” $\endgroup$
    – Dubukay
    May 9, 2019 at 15:17

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .