# An approximate, piecewise continuous analytical model for Earth's radial density profile useful for numerical integration?

An intriguing question in Space SE Optimal depth for underground flyby? asks about a theoretical spacecraft trajectory for a "slingshot maneuver" also known as a gravitational 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

update: This answer to Do there exist reasonable numbers on the density/radius/mass of the various geological layers of the Earth? 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.

• help with additional appropriate tags welcome!
– uhoh
Sep 23, 2022 at 1:31
• 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.) Sep 23, 2022 at 6:45
• @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
Sep 23, 2022 at 6:58
– uhoh
Sep 23, 2022 at 6:58
• @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. Sep 23, 2022 at 7:03

From Table I in the admittedly dated Preliminary Reference Earth Model,

$$\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.

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()

• 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. Sep 23, 2022 at 8:11
• Beautiful! I'll give it a test drive soon.
– uhoh
Sep 23, 2022 at 8:48

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/