Accurate Equation for Earth's Gravitational Binding Energy

This is a relatively important question for anyone who can answer it. I am trying to find the equation that accurately solves for Earth's Gravitational Binding Energy. The information below is from the wikipedia page:

Assuming that the Earth is a uniform sphere (which is not correct, but is close enough to get an order-of-magnitude estimate) with M = 5.97 x 10^24 kg and r = 6.37 x 10^6 m, U is 2.24 x 10^32 J. This is roughly equal to one week of the Sun's total energy output. It is 37.5 MJ/kg, 60% of the absolute value of the potential energy per kilogram at the surface.

The actual depth-dependence of density, inferred from seismic travel times (see Adams–Williamson equation), is given in the Preliminary Reference Earth Model (PREM). Using this, the real gravitational binding energy of Earth can be calculated numerically as U = 2.487 x 10^32 J.

So what I wish to not for the latter results (2.487e+32 J) is what is the actual equation is used to get this result.

And I do not mean the standard GBE equation which gave the other result above.

• In all likelihood they didn't use a closed-form equation, but instead numerically integrated a partial solution on the basis of the density information. This might have been a spherical calculation (in which case the thing the did the numeric integration on what what you called "the standard GBE equation"). If I was setting a challenge for a graduate mechanics class I would ask them to work it as a multi-pole problem as well, but that would probably add more uncertainty to the problem because it would require an assumption that the mutil-pole ratios were constant as a function of depth. – dmckee Nov 21 '18 at 15:56
• I found this which would likely get the answer: 1/(5 - c)*GBE , however I do not know the name of this equation, nor what the constant (c) value is – C. Jordan Nov 21 '18 at 16:09

To calculate the gravitational energy should use two integrals

$$M(r)=4\pi\int_0^r\rho(r')r'^2dr'$$

$$U=4\pi G\int_0^{R_E}\rho (r)M(r)rdr$$

$$R_E=6.371×10^6 m$$

We use the density distribution obtained on the basis of PREM (see Earth's Gravitational Binding Energy) As a result, we find $$U=2.48848×10^{32} J$$ which is consistent with the value $$2.487×10^{32} J$$ shown on the wikipedia page and calculated using PREM.

Below is the Wolfram code debugged for Mathematica 11.3.

A = {{0, 0, 0, 0}, {1.2215*10^6, -2.1773*10^-10, 1.9110*10^-8,
1.3088*10^4}, {3.4800*10^6, -2.4123*10^-10, 1.3976*10^-4,
1.2346*10^4}, {3.6300*10^6, 0, -5.0007*10^-4,
7.3067*10^3}, {5.7010*10^6, -3.0922*10^-11, -2.4441*10^-4,
6.7823*10^3}, {5.7710*10^6, 0, -2.3286*10^-4,
5.3197*10^3}, {5.9710*10^6, 0, -1.2603*10^-3,
1.1249*10^4}, {6.1510*10^6, 0, -5.9706*10^-4,
7.1083*10^3}, {6.3466*10^6, 0, 1.0869*10^-4,
2.6910*10^3}, {6.3560*10^6, 0, 0, 2.9000*10^3}, {6.3680*10^6, 0,
0, 2.6000*10^3}, {6.3710*10^6, 0, 0, 1.0200*10^3}};

rho[x_] :=
Piecewise[
Table[{A[[i, 2]]*x^2 + A[[i, 3]]*x + A[[i, 4]],
A[[i - 1, 1]] < x <= A[[i, 1]]}, {i, 2, Length[A]}]]
Plot[rho[x], {x, 0, 6.3710*10^6}, PlotRange -> All, PlotPoints -> 100,
AxesOrigin -> {0, 0},
AxesLabel -> {"r, m",
"\[Rho], kg/\!$$\*SuperscriptBox[\(m$$, $$3$$]\)"}]

G = 6.67408*10^\[Minus]11; R = 6.3710*10^6;
4*Pi*NIntegrate[
rho[x]*x^2, {x, 0,
R}](*a total Earth mass; the NASA figure of 5.9736\[Times]10^24 \
kg.*)

M =
Interpolation[
Table[{r, 4*Pi*NIntegrate[rho[x]*x^2, {x, 0, r}]}, {r, 1, R,
10^4}]];
Plot[M[x], {x, 0, R}, PlotRange -> All, AxesOrigin -> {0, 0},
AxesLabel -> {"r, m", "M, kg"}]

U =
4*Pi*G*NIntegrate[
rho[r]*r*M[r], {r, 0, R}](*Earth's Gravitational Binding Energy*)

(*Out[]= 2.48848*10^32*)

• Would there be anything else to know about any of the values in the integral equations you listed? – C. Jordan Nov 24 '18 at 15:21
• You can check that $M(R_E)=5.9727*10^{24} kg$ and compare with NASA data for the mass of the Earth $M_E=5.9736*10^{24} kg$. For the gravity constant, I used the value $G = 6.67408*10^{-11}$ from physics.nist.gov/cgi-bin/cuu/Value?bg . All calculations are done on Mathematica version 11.3 using the Piecewise[] function for $\rho (r)$. I can share the code. – Alex Trounev Nov 24 '18 at 16:19
• @C.Jordan I updated the message, added the code. – Alex Trounev Nov 27 '18 at 13:05
• Thank you Alex, and for Mathematica, is that a program you have to download? or is it free somewhere? – C. Jordan Nov 27 '18 at 15:51
• @C.Jordan You can use the trial version or install a cheap version for students wolfram.com/mathematica/?source=nav – Alex Trounev Nov 27 '18 at 16:00

I am trying to find the equation that accurately solves for Earth's Gravitational Binding Energy.

There isn't a single equation for this (and much of real science does not yield convenient single formulas as solutions). What the PREM produces for density is a set of piecewise approximate polynomial functions that model the theorized density that produces a reasonable match to measured data (like for example seismic data).

So working from this page by Dave Typinski you could develop an integral equation for the gravitational binding energy. I'm not going to actually do that myself, but the density polynomial functions are just quadratic, and you can apply this simply enough to the integral for gravitational binding energy :

$$U = 16\pi^2 G \int_0^R r \rho(r) \left[ \int_0^r \rho(r) r^2 dr \right] dr$$

This is tedious to do with the piecewise functions for density, but not difficult.

For a constant density $$\rho(r)=\rho_0$$ you can see this reduces to the familiar equation :

$$U = 16\pi^2G \rho_0^2 \int_0^R r \left[\frac 1 3 r^3\right] dr = \frac {16\pi^2G} {15} \rho_0^2 R^5 = \frac {3}{5} \frac {GM^2} R$$

So what I wish to not for the latter results ($$2.487\times 10^{32}\,J$$) is what is the actual equation is used to get this result.

I can't answer that definitively, but it would would most likely be either a numerical integration or something like I've outlined.