Modeling the Orbital Evolution of the Moon-Earth System

Revision as of 01:29, 30 April 2015 by Tenev (Talk | contribs)

Jump to: navigation, search



Dissipation of tidal energy in the Earth’s surface and interior as well as the pull of the Earth’s tidal bulge on the Moon causes the Moon’s orbit to recede away from the Earth at a rate that is currently about 3.8 cm/year. [1] This is consistent with model calculations such as Kaula, 1964 [2]. The rate of the lunar orbital recession has been used as an argument by Young-Earth creationists to place an upper limit on the age of the Earth-Moon system, which has brought about controversy due to doubts about the accuracy of the model used [3] [4]. Such model calculations are based on estimations of the specific energy dissipation function per unit volume Q, the bulk modulus k, and the rigidity modulus \mu throughout the Earth’s surface and interior, which do not lend themselves easily to direct measurements. In this paper we propose a method for computing these constitutive parameters by using ICME principles and performing simulations at multiple scales.


Figure 1. The force of attraction between the satellite m and the nearer tidal bulge A exceeds that between m and B; a component of the net torque retards the rotation of the planet M and accelerates the satellite in its orbit. [Included from Goldreich & Soter, 1965][5]

The Moon’s gravity causes the Earth to bulge approximately along the Earth-Moon orientation and to depress in the perpendicular plane. This applies not only to oceans, but also to the Earth’s crust and interior. However, because the Earth spins around its axis in one day, while it takes one month for the Moon to rotate in the same direction, the tidal bulge ends up slightly misaligned with the Earth-Moon direction as shown on Figure 1. Due to this asymmetry, the Earth would experience a body moment in the opposite direction from its rotation, i.e. in the direction of slowing down the Earth’s rotation, while at the same time the Moon would experience a body force component that would accelerate the Moon’s orbital velocity thus raising it up to a higher orbit. Effectively, a part of the Earth’s spin energy is being transferred to the Moon’s rotational energy. [6]

From the view point of someone on Earth near the equator, it would appear that the tidal bulge moves at a velocity of ≈1000 "miles/h" across the face of the Earth opposite to the direction of the Earth’s spin. The bending, compressing and tension of the crust and mantle would not be elastic and so some energy would be dissipated in the process. The bulge and energy dissipation are what ultimately determine the rate of the lunar orbital recession. These in turn depend on material specific parameters that can potentially be derived using ICME techniques.

Figure 2. Overview of ICME scales, simulations and bridges

The equation for the rate of lunar orbital recession as given by Kaula, 1964[2] is the following,

\dot{a} = - {{3 k_2 a_e^5 G m}\over{[G (M+m)]^{1/2} a^{11/2}}} sin{\epsilon_{2200}}

In the above equation, a is the semi-major axis of the Moon’s orbit; a_e is the mean radius of the Earth; M and m are the masses of the Earth and Moon, respectively, while G is the gravitational constant. The parameters \epsilon_{2200} and k_2 deserve more elaboration as explained below:

  • k_2 is a so called Love number (named after Augustus Edward Hough Love), which represents the correction that needs to be applied to the gravitational tidal potential as a result of the Earth’s own deformation in response to that potential. Thus, k_2 is also a kind of a measure of the rigidity of the Earth. For a perfectly rigid planet, k_2=0, but for Earth the actual value is k_2 \approx 0.3 [7]. Actually, this is an averaged value; in reality the value would vary for each point on the Earth’s globe and interior depending on its depth and location, such as land versus water or mantle versus core.
  • \epsilon_{2200} is the (0,0) component of the strain tensor in spherical coordinates. The subscripts (2,2) refer to the 2^nd degree and 2^nd order spherical harmonics (l=2,m=2). [8] That is, the \epsilon_{22}) strain tensor is one of the components of Fourier series expansion of the overall strain tensor, and that furthermore the (l=2,m=2) component has been determined to be the dominant one. [2]

In his paper Kaula, 1964[2], shows how one can calculate \epsilon_{2200} starting from the rigidity (a.k.a. shear) modulus \mu, bulk modulus k, of Earth, and the specific energy dissipation function Q. In the same paper, these values have been estimated and tabulated at different depths. The values for k_2 that Kaula uses are also based on estimates at different depths. Note that the specific energy dissipation function Q is defined as follows [5],

Q^{-1}={1\over {2 \pi E_0}} \oint {-{dE \over dt}}\,dt

where E_0 is the maximum energy stored in the tidal distortion and the integral over -dE/dt that is the rate of dissipation, is the energy lost during the complete cycle. Furthermore in Caputo, 1967 [9] the author offers an analytical method for computing Q based on the density \rho, shear modulus \mu, and Lamé’s first parameter \lambda, the latter being related to the shear modulus and the bulk modulus according to the following formula: \lambda=k-2\mu/3 [10] . Despite the proposed approximating formula for computing Q and k_2, the values of these parameters have been largely unknown for most bodies in the Solar system. Except for the Moon, for which the ratio k_2/Q=0.0011 is available and also for the Earth for which there are available estimates [2], for most other bodies the values for Q are arbitrarily chosen to be Q=100. [11][2]. The need to choose arbitrary values for Q shows the difficulty in obtaining these values just based on observation. Fundamentally, these are material dependent parameters that are internal state variables (ISVs).

We propose that the values of k_2 and Q be obtained through calculations at lower lengths scales. For example, we should be able to obtain Q from meter-length scale simulations using damage parameters obtained from lower scales for solids or using viscosity information for liquids. Likewise, we could calculate the various moduli from DFT simulations at the Angstrom scale. In turn, the damage parameters can be calculated from crystal plasticity simulations at the μm-scale and from dislocation density calculations at the Nm scale, while the viscosity information can be obtained from atomistics calculations.

To achieve a realistic model, the above simulations would be run for different parts of the Earth as follows:

  • Crust: assume composition predominantly of granite (for the land masses) and salt water (for oceans).
  • Mantle: assume composition of olivine.
  • Core: liquid iron and nickel
  • Inner core: solid iron.

Please refer to Figure 2 for an illustration of the pertinent length scale associated with each parameter.

Upscaling Methodology

Please refer to Figure 2 as the context for the following descriptions:

At the lowest length scales, we would run electronics principles calculations using density functional theory (DFT) techniques to obtain values for μ and k for the various materials and phases in which we are interested (salt water, granite, olivine, liquid/solid iron+nickel). This information would be passed to the macroscale (Bridge 4). In addition, the DFT simulation would also give us interfacial energy parameters which can be used in the higher length scale computation (Bridge 1) for obtaining values of viscosity (for liquids) or dislocation mobility parameters (for solids).

At the higher length scale, we would run atomistics simulations, e.g. Modified Embedded Atom Method (MEAM) to get kinematic parameters for atoms and molecules, which would feed mobility rules information into the dislocation dynamics simulation (Bridge 2). In the case of liquids, the atomistics simulation should also give us viscosity η (Bridge 5) to feed into the macroscale simulation.

The dislocation dynamics simulation should provide internal state variables (ISVs) characterizing dislocation density (Bridge 6) and hardening rules (Bridge 3).

Crystal plasticity finite element analysis (FEA) simulations would provide ISVs characterizing damage, i.e. voids, cracks nucleation, coalescence and propagation (Bridge 7). The information derived thus far should be sufficient to computationally determine the strain-stress response functions for each material of interest at the various pressures and temperatures throughout the globe.

Using the above mentioned ISVs, we could then run meter-scale FEA simulations to deduce the specific energy dissipation function Q as a function of the damage and dislocation density (Bridge 8).

Finally, armed with Q,\mu and k and the stress-strain response functions, we would be ready to run Terra 3D simulations at the scale of the entire globe to obtain the Love number k_2 at various depths and locations on the globe and, using information about the tidal body forces, to also compute the strain tensor \epsilon_{22ij} thus building up to the orbital evolution equation that was listed in the downscaling section above.


Once we have plumbed through all scale layers in this manner, we would be able to compute present-day values for the rate of orbital regression of the Moon and evaluate the uncertainty. These can be compared to actual values determined by laser measurements of the distance between the Earth and the Moon, such as described in Wikipedia, 2014 [12]. Furthermore the computations for Q and k_2 for the Earth at various depths as well as their respective uncertainties can be compared to current estimates such as the ones listed in Kaula, 1964 [2]


Figure 3. Model prediction of Earth-Moon distance in past years

Once we have constructed and validated the model in this way, we can use it to understand the behavior of the Earth-Moon system in the past. What should we expect? First, we could do a rough calculation based on the formula cited from Kaula, 1964[2] earlier. We can rewrite it in a simple form like this:

\dot{\alpha} \propto \alpha^{-11/2}

Knowing the current values for \dot{\alpha} = 3.8 \times 10^{-2} m/year and \alpha=3.84 \times 10^8 m we can compute the constant of proportionality to be:  6.25 10^45 m^{13/2}/year . From there we can run the calculation backwards and compute the Earth-Moon distance at different times in the past. The result is visualized on Figure 3.

As evidenced from the graph, the model appears "well behaved" up to about 1.5B years ago. For earlier times, the model predicts that due to the tidal effects examined here, the Moon would have to be regressing much faster from Earth and furthermore, it predicts that the whole process could not have been much older than 1.69B years at which point the Moon would have had to have been at a distance close to the geosynchronous orbit.

Young Earth creationists point to this result as an indication that the Moon could not have been older than 1.6B years in contrast to the popular science view that gives the age of 4.5B years. See for example Henry, 2006 or Bowden, 2000. Can the model be trusted? Could parameters that were negligibly small at larger Moon-Earth distances play a more significant role in the model at shorter distances and thus change the result? This is where a more thorough simulation with properly computed uncertainty bounds should help us arrive at a more definitive conclusion.


  1. Wikipedia (2015, 02 08). Lunar distance (astronomy). Retrieved 02 09, 2015, from Wikipedia:
  2. 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 Kaula, W.M. (1964, November). Tidal Dissipation by Solid Friction and the Resulting Orbital Evolution. Reviews of Geophysics, 2(4), 661-685
  3. Thompson, T. (1999, December). The recession of the Moon and the Age of the Earth-Moon System. Retrieved 02 09, 2015, from The TalkOrigins Archive:
  4. Bowden, M. (2000, 02-27). The Moon is Still Young - Rebutting Tim Thompson's "The Recession of the Moon" at Talk.Origins. Retrieved 02 09, 2015, from The True.Origin Archive - Exposing the Myth of Evolution:
  5. 5.0 5.1 Goldreich, P., & Soter, S. (1965, November 17). Q in the Solar System. Icarus, 375-389. Henry, J. (2006, August). The moon’s recession and age. Journal of Creation, 65-70.
  6. Pogge, R. (2007, October 14). Astronomy 161, Lecture 20:Tides. Retrieved 02 09, 2015, from Ohio State:
  7. Cite error: Invalid <ref> tag; no text was provided for refs named Poulsen09
  8. Wikipedia (2014, 11 28). Spherical Harmonics. Retrieved 02 09, 2015, from Wikipedia:
  9. Caputo, M. (1967). Linear Models of Dissipation whose Q is almost Frequency Independent-II. Geophysics Journal International, 529-539.
  10. Wikipedia (2015, January 15). Bulk Modulus. Retrieved 02 10, 2015, from Wikipedia:
  11. Gladman, B., Quinn, D., Nicholson, P., & Rand, R. (1996). Synchornous Locking of Tidally Evolving Satellites. Icarus, 166-192.
  12. Wikipedia (2014, December 31). Lunar Laser Ranging experiment. Retrieved 02 10, 2015, from Wikipedia:
Personal tools

Material Models