<p>From Table I in the admittedly dated <a href="https://doi.org/10.1016/0031-9201(81)90046-7" rel="noreferrer">Preliminary Reference Earth Model</a>,</p>
<p><span class="math-container">$$
\rho =
\begin{cases}
13.0885 - 8.8381 x^2 &amp; \quad \phantom{000}0\phantom{.0} \le r &lt; 1221.5 \\
12.5815 - 1.2638 x - 3.6426 x^2 - 5.5281 x^3 &amp; \quad 1221.5 \le r &lt; 3480.0 \\
\phantom{0}7.9565 - 6.4761 x + 5.5283 x^2 - 3.0807 x^3 &amp; \quad 3480.0 \le r &lt; 5701.0 \\
\phantom{0}5.3197 - 1.4836 x &amp; \quad 5701.0 \le r &lt; 5771.0 \\
11.2494 - 8.0298 x &amp; \quad 5771.0 \le r &lt; 5971.0 \\
\phantom{0}7.1089 - 3.8045 x &amp; \quad 5971.0 \le r &lt; 6151.0 \\
\phantom{0}2.6910 + 0.6924 x &amp; \quad 6151.0 \le r &lt; 6346.6 \\
\phantom{0}2.900 &amp; \quad  6346.6 \le r &lt; 6356.0 \\
\phantom{0}2.600 &amp; \quad  6356.0 \le r &lt; 6368.0 \\
\phantom{0}1.020 &amp; \quad 6291.0 \le r \le 6371.0 \\
\end{cases}
$$</span>
where</p>
<ul>
<li><span class="math-container">$\rho$</span> is the density in grams per cubic centimeter,</li>
<li><span class="math-container">$r$</span> is distance from the center of the Earth to the point of interest in kilometers, and</li>
<li><span class="math-container">$x$</span> is the normalized radius: <span class="math-container">$x\equiv\frac r {6371}$</span>.</li>
</ul>
<p>As each element is a polynomial, this is a piecewise continuous model of density inside the Earth from the center to the surface. There are however discontinuities at the transitions.</p>
<p>A plot of PREM density as a function of radial distance is portrayed below.
<a href="https://i.sstatic.net/Kq1aV.png" rel="noreferrer"><img src="https://i.sstatic.net/Kq1aV.png" alt="PREM density model, density in g/cc as a function of radial distance in km" /></a></p>
<p>I generated this using the following python script:</p>
<pre><code>from matplotlib import pyplot as plt

def density(r) :
    x = r / 6371.0
    if r &lt; 1221.5 :
        rho = 13.0885 - x*x*8.8381
    elif r &lt; 3480.0 :
        rho = 12.5815 - x*(1.2638 + x*(3.6426 + x*5.5281))
    elif r &lt; 5701.0 :
        rho = 7.9565 - x*(6.4761 - x*(5.5283 - x*3.0807))
    elif r &lt; 5771.0 :
        rho = 11.2494 - x*8.0298
    elif r &lt; 6151.0 :
        rho = 7.1089 - x*3.8045
    elif r &lt; 6346.6 :
        rho = 2.6910 + x*0.6924
    elif r &lt; 6356.0 :
        rho = 2.900
    elif r &lt; 6368.0 :
        rho = 2.600
    else :
        rho = 1.020
    return rho

def main() :
    r = [float(x) for x in range(6372)]
    rho = [density(x) for x in r]
    plt.plot(r, rho)
    plt.ylabel('rho (g/cc)')
    plt.xlabel('r (km)')
    plt.title('PREM density')
    plt.grid()
    plt.savefig('prem_dens.png')
    plt.show()

if __name__ == &quot;__main__&quot; :
    main()
</code></pre>

CC BY-SA 4.0
2022-09-23T08:11:54.350Z
I tripled checked the plus sign in $2.6910+0.6924x$ for the low velocity zone (LVZ) and the rigid lithosphere (LID). It's definitely a plus sign rather than a minus sign. Per this model, density increases with increasing radial distance (decreases with increasing depth) in those regions.
David Hammen
CC BY-SA 4.0
2022-09-23T08:48:59.510Z
Beautiful! I'll give it a test drive soon.
uhoh
CC BY-SA 4.0
2022-09-23T07:52:07.757Z
2022-09-23T09:20:53.960Z
2022-09-23T09:20:53.960Z
David Hammen
<p>If you want access to the raw data, the Preliminary Reference Earth Model (and many other more refined models) are available on the IRIS website:</p>
<p><a href="https://ds.iris.edu/ds/products/emc-prem/" rel="nofollow noreferrer">https://ds.iris.edu/ds/products/emc-prem/</a></p>
<p>For example you can download the file PREM_1s_IDV.csv</p>
<p>Each line represents a depth and a density value at this depth (columns &quot;depth and &quot;density&quot;). From this you can construct a piecewise discontinuous linear model, or directly use them for a numerical integration.</p>

CC BY-SA 4.0
2022-09-26T15:53:28.963Z
2022-09-26T15:53:28.963Z
Navid
<p>An intriguing question in Space SE <a href="https://space.stackexchange.com/q/60417/12102">Optimal depth for underground flyby?</a> asks about a theoretical spacecraft trajectory for a &quot;slingshot maneuver&quot; also known as a <a href="https://en.wikipedia.org/wiki/Gravity_assist" rel="noreferrer">gravitational assist</a> or &quot;flyby&quot; with the relaxed condition that the trajectory can go through the earth (presumably through a carefully constructed tunnel who's path is adjusted to account for Earth's rotation during the f̶l̶y̶b̶y̶ &quot;flythrough&quot;.</p>
<p>There are answers there but currently the most rigorous numerical trajectories still assume a uniform density sphere. The answer will be much different if it can include the dramatic <em>factor of 5</em> radial variation of Earth's density from say 2.5 (ignoring water) to over 13 g/cm<sup>3</sup>.</p>
<p>I found the plot below in Wikimedia which certainly looks like it's piecewise polynomial and perhaps even piecewise parabolic.</p>
<p>Sadly, it references only <a href="https://www.gps.caltech.edu/uploads/File/People/dla/DLApepi81.pdf" rel="noreferrer">https://www.gps.caltech.edu/uploads/File/People/dla/DLApepi81.pdf</a> as the source of the data and that link is of course broken (ones like this usually do) rather than a proper citation (name, year, journal, page number) or a doi link.</p>
<p>So I'd like to ask the following.</p>
<p><strong>Question:</strong> Is there an approximate, piecewise continuous analytical model for Earth's radial density profile useful for numerical integration that can be used as a properly cite-able source for a Stack Exchange answer?</p>
<p>One could reverse-engineer this one, but it would be better to get at the source.</p>
<hr />
<p><a href="https://i.sstatic.net/35yyl.jpg" rel="noreferrer"><img src="https://i.sstatic.net/35yyl.jpg" alt="Earth's radial density distribution according to the Preliminary Reference Earth Model (PREM). 29 December 2010" /></a></p>
<blockquote>
<p>Earth's radial density distribution according to the Preliminary Reference Earth Model (PREM). 29 December 2010</p>
</blockquote>
<p>Source: <a href="https://commons.wikimedia.org/wiki/File:RadialDensityPREM.jpg" rel="noreferrer">Wikimedia's File: RadialDensityPREM.jpg</a></p>
<p><strong>update:</strong> <a href="https://earthscience.stackexchange.com/a/15644/6031">This answer to <em>Do there exist reasonable numbers on the density/radius/mass of the various geological layers of the Earth?</em></a> provides one such model for Earth's core with two polynomials, but these end at a radius of 3480 km with nothing outside of that.</p>

CC BY-SA 4.0
2022-09-23T01:31:34.080Z
help with additional appropriate tags welcome!
uhoh
CC BY-SA 4.0
2022-09-23T06:45:17.367Z
Frame challenge: (1) Those discontinuities are real. Why do want to wash them away? (2) Numerical integration across a discontinuity is not hard. (Differentiation is a challenge, but not integration.)
David Hammen
CC BY-SA 4.0
2022-09-23T06:58:01.770Z
@DavidHammen I've asked for "An approximate, continuous piecewise analytical model..." it doesn't have to be continuous everywhere, just have a manageable number of regions. Seems like five or six would be fine. Don't let Fred's mischaracterization fo my question fool you.
uhoh
CC BY-SA 4.0
2022-09-23T06:58:24.847Z
@Fred please read carefully before writing misdirections.
uhoh
CC BY-SA 4.0
2022-09-23T07:03:06.183Z
@uhoh I didn't base my comment on Fred's comment. I based it on the title to the question. If you want a somewhat realistic answer to that non-intriguing (to me) question on SE, use the PREM. There are ways to get your hands on that paper (illegal in some countries). That paper gives piecewise continuous polynomial expressions for various things, including gravitational acceleration, inside the Earth. But there are discontinuities because the discontinuities are real.
David Hammen
CC BY-SA 4.0
2022-09-23T07:05:39.290Z
@DavidHammen what would you change about the title to make it clearer to you, perhaps "piecewise continuous" instead of "continuous piecewise"?
uhoh
CC BY-SA 4.0
2022-09-23T07:13:34.587Z
I suspect you know how to get that paper without paying the paywall fee. if you don't, google Alexandra Elbakyan.
David Hammen
CC BY-SA 4.0
2022-09-23T07:17:19.933Z
Changing the title to 'piecewise continuous' would help. The order of the words matter in this case.
David Hammen
CC BY-SA 4.0
2022-09-23T08:46:39.373Z
@DavidHammen done. I'll stop by campus and see if I can locate the paper. *Thanks!*
uhoh
<open-data><structural-geology>
An approximate, piecewise continuous analytical model for Earth's radial density profile useful for numerical integration?
uhoh
# An approximate, piecewise continuous analytical model for Earth's radial density profile useful for numerical integration?

(Asked by: username_0 on 2022-09-23)

An intriguing question in Space SE [Optimal depth for underground flyby?](https://space.stackexchange.com/q/60417/12102) asks about a theoretical spacecraft trajectory for a "slingshot maneuver" also known as a [gravitational assist](https://en.wikipedia.org/wiki/Gravity_assist) or "flyby" with the relaxed condition that the trajectory can go through the earth (presumably through a carefully constructed tunnel who's path is adjusted to account for Earth's rotation during the f̶l̶y̶b̶y̶ "flythrough".

There are answers there but currently the most rigorous numerical trajectories still assume a uniform density sphere. The answer will be much different if it can include the dramatic *factor of 5* radial variation of Earth's density from say 2.5 (ignoring water) to over 13 g/cm3.

I found the plot below in Wikimedia which certainly looks like it's piecewise polynomial and perhaps even piecewise parabolic.

Sadly, it references only <https://www.gps.caltech.edu/uploads/File/People/dla/DLApepi81.pdf> as the source of the data and that link is of course broken (ones like this usually do) rather than a proper citation (name, year, journal, page number) or a doi link.

So I'd like to ask the following.

**Question:** Is there an approximate, piecewise continuous analytical model for Earth's radial density profile useful for numerical integration that can be used as a properly cite-able source for a Stack Exchange answer?

One could reverse-engineer this one, but it would be better to get at the source.

---

[![Earth's radial density distribution according to the Preliminary Reference Earth Model (PREM). 29 December 2010](https://i.sstatic.net/35yyl.jpg)](https://i.sstatic.net/35yyl.jpg)

> Earth's radial density distribution according to the Preliminary Reference Earth Model (PREM). 29 December 2010

Source: [Wikimedia's File: RadialDensityPREM.jpg](https://commons.wikimedia.org/wiki/File:RadialDensityPREM.jpg)

**update:** [This answer to *Do there exist reasonable numbers on the density/radius/mass of the various geological layers of the Earth?*](https://earthscience.stackexchange.com/a/15644/6031) provides one such model for Earth's core with two polynomials, but these end at a radius of 3480 km with nothing outside of that.

**Tags:** open-data, structural-geology

> **Comments:**
> **username_0**: help with additional appropriate tags welcome!  
> **username_1**: Frame challenge: (1) Those discontinuities are real. Why do want to wash them away? (2) Numerical integration across a discontinuity is not hard. (Differentiation is a challenge, but not integration.)  
> **username_0**: @username_1 I've asked for "An approximate, continuous piecewise analytical model..." it doesn't have to be continuous everywhere, just have a manageable number of regions. Seems like five or six would be fine. Don't let Fred's mischaracterization fo my question fool you.  
> **username_0**: @Fred please read carefully before writing misdirections.  
> **username_1**: @username_0 I didn't base my comment on Fred's comment. I based it on the title to the question. If you want a somewhat realistic answer to that non-intriguing (to me) question on SE, use the PREM. There are ways to get your hands on that paper (illegal in some countries). That paper gives piecewise continuous polynomial expressions for various things, including gravitational acceleration, inside the Earth. But there are discontinuities because the discontinuities are real.  
> **username_0**: @username_1 what would you change about the title to make it clearer to you, perhaps "piecewise continuous" instead of "continuous piecewise"?  
> **username_1**: I suspect you know how to get that paper without paying the paywall fee. if you don't, google Alexandra Elbakyan.  
> **username_1**: Changing the title to 'piecewise continuous' would help. The order of the words matter in this case.  
> **username_0**: @username_1 done. I'll stop by campus and see if I can locate the paper. *Thanks!*  

---

## 2 answers

### Answer by username_1 on 2022-09-23. Score: 7

From Table I in the admittedly dated [Preliminary Reference Earth Model](https://doi.org/10.1016/0031-9201(81)90046-7),

$$
\rho =
\begin{cases}
13.0885 - 8.8381 x^2 & \quad \phantom{000}0\phantom{.0} \le r < 1221.5 \\
12.5815 - 1.2638 x - 3.6426 x^2 - 5.5281 x^3 & \quad 1221.5 \le r < 3480.0 \\
\phantom{0}7.9565 - 6.4761 x + 5.5283 x^2 - 3.0807 x^3 & \quad 3480.0 \le r < 5701.0 \\
\phantom{0}5.3197 - 1.4836 x & \quad 5701.0 \le r < 5771.0 \\
11.2494 - 8.0298 x & \quad 5771.0 \le r < 5971.0 \\
\phantom{0}7.1089 - 3.8045 x & \quad 5971.0 \le r < 6151.0 \\
\phantom{0}2.6910 + 0.6924 x & \quad 6151.0 \le r < 6346.6 \\
\phantom{0}2.900 & \quad 6346.6 \le r < 6356.0 \\
\phantom{0}2.600 & \quad 6356.0 \le r < 6368.0 \\
\phantom{0}1.020 & \quad 6291.0 \le r \le 6371.0 \\
\end{cases}
$$
where

* $\rho$ is the density in grams per cubic centimeter,
* $r$ is distance from the center of the Earth to the point of interest in kilometers, and
* $x$ is the normalized radius: $x\equiv\frac r {6371}$.

As each element is a polynomial, this is a piecewise continuous model of density inside the Earth from the center to the surface. There are however discontinuities at the transitions.

A plot of PREM density as a function of radial distance is portrayed below.
[![PREM density model, density in g/cc as a function of radial distance in km](https://i.sstatic.net/Kq1aV.png)](https://i.sstatic.net/Kq1aV.png)

I generated this using the following python script:

```
from matplotlib import pyplot as plt

def density(r) :
    x = r / 6371.0
    if r < 1221.5 :
        rho = 13.0885 - x*x*8.8381
    elif r < 3480.0 :
        rho = 12.5815 - x*(1.2638 + x*(3.6426 + x*5.5281))
    elif r < 5701.0 :
        rho = 7.9565 - x*(6.4761 - x*(5.5283 - x*3.0807))
    elif r < 5771.0 :
        rho = 11.2494 - x*8.0298
    elif r < 6151.0 :
        rho = 7.1089 - x*3.8045
    elif r < 6346.6 :
        rho = 2.6910 + x*0.6924
    elif r < 6356.0 :
        rho = 2.900
    elif r < 6368.0 :
        rho = 2.600
    else :
        rho = 1.020
    return rho

def main() :
    r = [float(x) for x in range(6372)]
    rho = [density(x) for x in r]
    plt.plot(r, rho)
    plt.ylabel('rho (g/cc)')
    plt.xlabel('r (km)')
    plt.title('PREM density')
    plt.grid()
    plt.savefig('prem_dens.png')
    plt.show()

if __name__ == "__main__" :
    main()

```

> **Comments:**
> **username_1**: I tripled checked the plus sign in $2.6910+0.6924x$ for the low velocity zone (LVZ) and the rigid lithosphere (LID). It's definitely a plus sign rather than a minus sign. Per this model, density increases with increasing radial distance (decreases with increasing depth) in those regions.  
> **username_0**: Beautiful! I'll give it a test drive soon.  

---

### Answer by username_2 on 2022-09-26. Score: 2

If you want access to the raw data, the Preliminary Reference Earth Model (and many other more refined models) are available on the IRIS website:

<https://ds.iris.edu/ds/products/emc-prem/>

For example you can download the file PREM\_1s\_IDV.csv

Each line represents a depth and a density value at this depth (columns "depth and "density"). From this you can construct a piecewise discontinuous linear model, or directly use them for a numerical integration.


---

