We can estimate how density varies inside the Earth by using seismic data, particularly the speeds of P-waves and S-waves. These waves travel through the Earth’s interior during earthquakes and carry information about its internal structure.

The velocities of seismic waves depend on the elastic properties and density of the material.

• P-wave (primary, compressional) speed:

$$ v_p^2 = \frac{K + \frac{4}{3} \mu}{\rho} $$

• S-wave (secondary, shear) speed:

$$ v_s^2 = \frac{\mu}{\rho} $$

Here:

• $K$ = bulk modulus (resistance to compression)
• $\mu$ = shear modulus (resistance to shear)
• $\rho$ = density

By eliminating $\mu$, we can express the bulk modulus per unit density as:

$$ \frac{K}{\rho} = v_p^2 - \frac{4}{3} v_s^2 $$

The pressure at a given depth in the Earth must balance the weight of overlying layers. The pressure difference between the top and bottom of a spherical shell is given by:

$$ \frac{\Delta P}{\Delta r} \simeq -\frac{G \rho M(<r)}{r^2} $$

where:

• $G$ = gravitational constant
• $M(<r)$ = mass enclosed within radius $r$

The bulk modulus is also defined as:

$$ K = \frac{\Delta P}{\Delta \rho} $$

Combining the two equations, we get the density gradient equation:

$$ \Delta \rho = - \frac{\rho^2 G M(<r)}{K r^2} \Delta r $$

This tells us how density changes with depth, provided we know seismic wave speeds (to get $K/\rho$) and mass distribution.

To compute the density profile from the surface to the center of the Earth:

1. Step 1: Use P- and S-wave travel time data to calculate $K/\rho$ at each depth.
2. Step 2: Assume outermost density (surface rock) $\rho \approx 3000$ kg/m³.
3. Step 3: Estimate total mass $M = 5.972 \times 10^{24}$ kg.
4. Step 4: Discretize the Earth into thin shells (e.g., 100 layers).
5. Step 5: Use the density gradient equation to compute $\Delta \rho$ for each shell.
6. Step 6: Recompute the mass enclosed within each radius.
7. Step 7: Repeat inward to the core.

Here’s a Python implementation of the method:

```python import numpy as np import matplotlib.pyplot as plt

# Constants G = 6.67430e-11 R_earth = 6.371e6 rho_surface = 3000 M_earth = 5.972e24 N = 100 dr = R_earth / N

# Seismic speeds vp = 8000 # m/s vs = 4500 # m/s K_over_rho = vp2 - (4/3) * vs2

# Initialize r = np.linspace(R_earth, 0, N) rho = np.zeros(N) mass = np.zeros(N) rho[0] = rho_surface mass[0] = M_earth

# Loop for i in range(1, N):

  mass[i] = mass[i-1] - 4*np.pi*r[i-1]**2 * dr * rho[i-1]
  delta_rho = - (rho[i-1]**2 * G * mass[i]) / (K_over_rho * r[i-1]**2) * dr
  rho[i] = rho[i-1] + delta_rho

# Plot plt.plot(r / 1000, rho) plt.xlabel(“Radius (km)”) plt.ylabel(“Density (kg/m³)”) plt.title(“Earth’s Density Profile (Williams–Adams Algorithm)”) plt.gca().invert_xaxis() plt.show()