View source for Modeling the Orbital Evolution of the Moon-Earth System

; Abstract

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.
<ref name=wiki_a> Wikipedia (2015, 02 08). Lunar distance (astronomy). Retrieved 02 09, 2015, from Wikipedia: http://en.wikipedia.org/wiki/Lunar_distance_(astronomy)</ref> 
This is consistent with model calculations such as Kaula, 1964
<ref name=Kaula64> Kaula, W.M. (1964, November). Tidal Dissipation by Solid Friction and the Resulting Orbital Evolution. ''Reviews of Geophysics'', 2(4), 661-685</ref>. 
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 
<ref name=Thompson99> 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: http://www.talkorigins.org/faqs/moonrec.html </ref>
<ref name=Bowden00> 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: http://www.trueorigin.org/moonmb.asp </ref>. 
Such model calculations are based on estimations of the specific energy dissipation function per unit volume <math>Q</math>, the bulk modulus <math>k</math>, and the rigidity modulus <math>\mu</math> 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.

== Background ==

[[Image:Earth Tidal Buldge.jpg|400px|thumb| 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]<ref name=Goldreich65> 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. </ref>]]

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. <ref name=Pogge07> Pogge, R. (2007, October 14). Astronomy 161, Lecture 20:Tides. Retrieved 02 09, 2015, from Ohio State: http://www.astronomy.ohio-state.edu/~pogge/Ast161/Unit4/tides.html </ref>

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.

[[Image:Earth-Moon_ICME_Scales.jpg|400px|thumb| Figure 2. Overview of ICME scales, simulations and bridges]]

The equation for the rate of lunar orbital recession as given by Kaula, 1964<ref name=Kaula64/> is the following,

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

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

* <math>k_2</math> 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, <math>k_2</math> is also a kind of a measure of the rigidity of the Earth. For a perfectly rigid planet, <math>k_2=0</math>, but for Earth the actual value is <math>k_2 \approx 0.3</math> <ref name=Poulsen09>Poulsen, S. K. (2009). ''Tidal Deformatoin of the Solid Earth - A Finite Difference Discretization.'' Faculty of Science. Niels Bohr Institute, University of Copenhagen.</ref>. 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.
	
* <math>\epsilon_{2200}</math> 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). <ref name=wiki_b> Wikipedia (2014, 11 28). Spherical Harmonics. Retrieved 02 09, 2015, from Wikipedia: http://en.wikipedia.org/wiki/Spherical_harmonics </ref> That is, the <math>\epsilon_{22})</math> 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. <ref name=Kaula64/>

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

:<math>
Q^{-1}={1\over {2 \pi E_0}} \oint {-{dE \over dt}}\,dt
</math>

where <math>E_0</math> is the maximum energy stored in the tidal distortion and the integral over <math>-dE/dt</math> that is the rate of dissipation, is the energy lost during the complete cycle. Furthermore in Caputo, 1967 <ref name=Caputo67> Caputo, M. (1967). Linear Models of Dissipation whose Q is almost Frequency Independent-II. Geophysics Journal International, 529-539. </ref>
the author offers an analytical method for computing <math>Q</math> based on the density <math>\rho</math>, shear modulus <math>\mu</math>, and Lamé’s first parameter <math>\lambda</math>, the latter being related to the shear modulus and the bulk modulus according to the following formula: <math>\lambda=k-2\mu/3</math> <ref name=wiki_c>Wikipedia (2015, January 15). Bulk Modulus. Retrieved 02 10, 2015, from Wikipedia: http://en.wikipedia.org/wiki/Bulk_modulus </ref>
. Despite the proposed approximating formula for computing <math>Q</math> and <math>k_2</math>, the values of these parameters have been largely unknown for most bodies in the Solar system. Except for the Moon, for which the ratio <math>k_2/Q=0.0011</math> is available and also for the Earth for which there are available estimates <ref name=Kaula64/>, for most other bodies the values for <math>Q</math> are arbitrarily chosen to be <math>Q=100</math>. <ref name=Gladman96> Gladman, B., Quinn, D., Nicholson, P., & Rand, R. (1996). Synchornous Locking of Tidally Evolving Satellites. Icarus, 166-192. </ref><ref name=Kaula64/>. The need to choose arbitrary values for <math>Q</math> 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 <math>k_2</math> and <math>Q</math> be obtained through calculations at lower lengths scales. For example, we should be able to obtain <math>Q</math> 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 <math>Q,\mu</math> and <math>k</math> 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 <math>k_2</math> at various depths and locations on the globe and, using information about the tidal body forces, to also compute the strain tensor <math>\epsilon_{22ij}</math> thus building up to the orbital evolution equation that was listed in the downscaling section above.

== Validation ==

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
<ref name=wiki_d>Wikipedia (2014, December 31). Lunar Laser Ranging experiment. Retrieved 02 10, 2015, from Wikipedia: http://en.wikipedia.org/wiki/Lunar_Laser_Ranging_experiment </ref>.
Furthermore the computations for <math>Q</math> and <math>k_2</math> 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 <ref name=Kaula64/>

== Application ==
[[Image:Earth-Moon_Distance.jpg|thumb|400px|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<ref name=Kaula64/> earlier. We can rewrite it in a simple form like this:

:<math>
\dot{\alpha} \propto \alpha^{-11/2}
</math>

Knowing the current values for <math>\dot{\alpha} = 3.8 \times 10^{-2} m/year</math> and <math>\alpha=3.84 \times 10^8 m </math>  we can compute the constant of proportionality to be: <math> 6.25 10^45 m^{13/2}/year </math> . 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.

== Future Work ==
The proposal above is limited to the calculation of the orbital evolution of the Earth-Moon system, but the modeling work that will be involved would enable related investigations, such as relating to the effects of the tidal energy dissipation on the Earth itself. The latter too lends itself to be investigated by the ICME techniques and is in fact quite significant. For example Denis, et al, 2002 <ref name=Denis02>
Denis, C., Schreider, A., Varga, P., & Zavoti, J. (2002). Despinning of the earth rotation in the geological past and. ''Journal of Geodynamics'', 667-685.
</ref>
quote Maslov, 1991 <ref name=Maslov91>
Maslov, L. (1991). Geodynamics of the Pacific segment of the Earth. ''Nauka (Science) in Russian''.
</ref>
and provide the following comparison:

{|class="wikitable"
|-
!Type of Energy	
! Amount [J/a]
|-
| Solar energy received by the Earth	
| <math>2.1 \times 10^{24} </math>
|-
| Loss of energy by heat flow	
| <math>1.0 \times 10^{19} </math>
|-
| Dissipation of energy by tidal friction	
| <math>1.6 \times 10^{19} </math>
|-
| Energy released by tectonic activity	
| <math>1.3 \times 10^{19} </math>
|-
| Energy released by seismic waves	
| <math>1.0 \times 10^{18} </math> 
|}

The above shows that the tidal energy dissipation per year is comparable to the energies of other Earth-shaping phenomena and so it is likely to have significant effects especially in geological areas of high internal stresses due to their specific geometry or formation history. For example Su, et al, 2012 <ref name=Su12>
Su, Y., Fu, H., & Hu, H. (2012). Study on Correlation of Tidal Forces with Global Great. ''International Journal of Geosciences'', 373-378.
</ref>
point out the long history of study of the correlation between tidal forces and seismic activity and conclude that “the lunar node tide is possibly one of the important astronomical factors to the seismic activities [in the principal seismic belts and regions of the world]”. Presumably, through FEA modeling of the Earth’s globe and the tidal effects, we should be able to make a better qualification and quantification of this assertion.
	
As part of such work, we could also consider the possibility for tidal energy to accumulate in the form of useful potential energy that is not necessarily released at the end of a tidal cycle but rather builds up and later released abruptly thereby causing local or perhaps even global catastrophes. A similar idea had been suggested by Brown 2008
<ref name=Brown08>
Brown, W. (2008). ''In the Beginning: Compelling Evidence for Creation and the Flood (8th Edition)''. Center for Scientific Creation.
</ref>
in his Hydroplate Theory of The Flood. Although the theory in its entirety has been a subject of much controversy and debunking from both Young Earth creationists as well as secular scientists, the description of the tidal pumping action is an example of an effect that could lead to the accumulation of usable potential energy. If an analogous, albeit not identical, mechanism exists today on Earth, its discovery could lead to a greater understanding of the Earth’s geological past as well as perhaps better prediction of future cataclysms.

== References ==

<references />