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

Question

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/cm³.

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

Source: Wikimedia's 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? 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.

Answer 1

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

Answer 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.