Abstract

Paleoclimate data have yielded variations with periods of ~23, ~40, and ~100 ky. Thermodynamic changes resulting from orbital eccentricity, obliquity, and precession have been ascribed as the cause of the variations although processes within the oceans and atmosphere may have too short memory to explain such variations. In this work, the dynamics of Sun-Moon gravitation (SMG) were explored for a rotating Earth and were determined to have a long memory in magma, a mostly ignored geophysical fluid with a mass ~3,400 times that of the atmosphere plus the oceans. Using the basic motion and gravitation (including obliquity) of the Sun and the Moon, we determined that SMG-induced magma motion could produce paleoclimatic variations with multiple periods (e.g., ~23, ~40, ~80, and ~100 ky), with considerable power for Earth’s heat. Such “reproducible” power could possibly maintain an energetic Earth against collapse, radioactivity, and cooling.

1. Introduction

Orbital eccentricity, obliquity, and precession, with periodicities of ~100, ~41, and ~23 ky, respectively, are three important drivers of paleoclimatic variations [1–4]. Each type of orbital forcing thermodynamically impacts climate, by altering either the total amount or the distribution of insolation received on Earth. For example, latitude-averaged insolation is reduced during Northern Hemisphere summer, leading to the reduction of summer snow-melt at high latitudes and over long time periods to glacial expansion due to lower obliquity [5]. Ice core records obtained from Antarctica indicate that orbitally driven Antarctic climate change lags behind Northern Hemisphere insolation variations by a few millennia [6, 7]. Insolation (i.e., Milankovitch Cycles) summarizes the thermodynamics of active orbital forcing. However, as indicated in Figure 1, we determined that insolation yielded a low correlation (~24.9%) to and a low contribution to temperature on paleoclimatic scale and yielded different anomalies and spectra from those of observed temperature and CO2 concentrations over the past 800 ky [6]. The observed temperature and CO2 concentrations over the past 800 ky highly correlated to each other, with an ~88.6% correlation with a 99% confidence. The spectra of both temperature and CO2 displayed periods of ~23 (22.30–23.73) ky, ~40 (38.79–41.27) ky, ~70 (67.88–72.22) ky, and ~100 (96.97–103.17) ky, while the spectra of insolation produced shorter periods of ~20–23 (19.39–23.73) ky and ~41 (39.75–42.30) ky. Period ranges within the parentheses above were from error estimations via a chi-square test [8, 9] for a 95% significance level. Also, temperature and CO2 yielded a much longer total phase during negative-anomaly phases (443 and 412 ky) than during positive-anomaly phases (357 and 389 ky), while insolation had a shorter total phase during negative-anomaly phases (391 ky) than during positive-anomaly phases (409 ky). Finally, temperature and CO2 yielded much smaller amplitudes during negative-anomaly phases (~5°C and ~52 ppm) than during positive-anomaly phases (~10°C and ~75 ppm), while insolation had approximately the same amplitude during negative- and positive-anomaly phases (~63 and ~66 Wm−2, resp.). Both temperature and CO2 yielded a “buffer” or a slower accumulation during negative anomaly phases, but a relatively quick release during positive anomaly phases, while insolation did not.

Figure 1: A time series with the mean deducted (left column, data interval is 0.1 ky) and its squared amplitude spectrum (right column) for temperature (Row a), CO2 (Row b), insolation (i.e., Milankovitch cycle, Row c), and dust (Row d) from the EPICA Ice Core during July at 65°N for the past 800,000 years [6, 10, 11], as well as the total PSMGIM (1 TW = 1012 W, Row e, using cf = 0.001 s−1, , and 4 μm/s; for the lower mantle and outer core, resp., see Section 2). Left column: the time series (red-green curve) was 100% correlated to its reconstruction (dashed blue curve) that was produced via the Discrete Fourier Transform (see Methods, Section 2.3), with a confidence of 99.99%, and was evaluated using the following: the positive/negative amplitude , the total phase length during positive-/negative-phases , the mean , and the correlation to temperature (CT). Right column: the squared amplitude spectrum (see Methods, Section 2.3, black curve) with the -axis representing the value and the -axis representing the period (ky). Error bars were plotted with a dashed blue (upper limit) and red (lower limit) line, estimated via a chi-square test [8, 9] for a 95% significance level.

Past studies on slow processes within the atmosphere and oceans have indicated climatic variations of shorter periods from ~1 to ~100 ky [12–16]. For example, the approximate 1.5 ky period in the north Atlantic Ocean climate system is caused by a combination of weak periodic forcing and “noise” from ice sheet related events [12]. The 1 ky Northern Hemisphere temperature and the 100 ky cycle for tropical sea surface temperature are associated with greenhouse gases [13, 14]. Regional climate shifts during the Pliocene epoch (~4.5 to 3.0 million years ago) are related to gradual global cooling [15]. Changes in sea level and temperature are also influenced by freshwater forcing via the phase change [16].

However, the slow processes within the atmosphere and oceans may only explain shorter variation probably because the atmosphere and oceans have memories shorter than those of paleoclimatic variations from the viewpoint of dynamics [17], and processes within the atmosphere and oceans may largely serve as the responsive portion of paleoclimatic variations. Slower driver(s) should exist that can influence CO2 and temperature on paleoclimate scales. Ice core data obtained from Antarctica over the past 800 ky provide detailed insights into the aerosol (dust) load of the atmosphere [10]. We determined that dust is highly correlated with temperature, with a correlation coefficient of −67.5% (confidence > 97%, via a double-sided -test). Dust lags temperature by approximately 500 years and has almost the same spectra as temperature, with periods of ~23–25 (22.30–25.79) ky, ~40 (38.79–41.27) ky, ~75 (72.73–77.38) ky, and ~100 (96.97–103.17) ky (from Figure 1(d), the period ranges within the parentheses above were estimated via a chi-square test [8, 9] for a 95% significance level). During the period of time from 800 ky (BP) to the preindustrial period, the dust load should largely result from natural processes associated with volcanism and, therefore, should be related to the activity of magma associated with variations in temperature and dust.

Using a theoretical model established for SMG-induced magma motion (see Section 2), we examined the Power of SMG-Induced Magma (PSMGIM) motion and determined that PSMGIM results in significant variations on multiple paleoclimatic scales for an 800 ky time window that began from the time period set in the model (instead of 800 ky BP). Our experimental solutions for the PSMGIM yielded spectra similar to those for temperature and CO2 for “buffer” and “quick release” mechanisms during negative and positive anomaly phases, respectively. The spectra of the PSMGIM displayed periods of ~40 (38.79–41.27) ky, ~80 (77.58–82.54) ky, and ~100 (96.97–103.17) ky. Period ranges within the parentheses above were estimated via a chi-square test [8, 9] for a 95% significance level. As shown in Figure 1(e), the total phase length was much longer during negative anomaly phases (417 ky) than during positive anomaly phases (383 ky), and the amplitude was much smaller during negative anomaly phases (1,206 TW, 1 TW = 1012 W) than during positive anomaly phases (2,063 TW). Shorter periods (e.g., of ~1-2 and ~22.3–23.7 ky) were also determined to appear when smaller diffusion coefficients were employed (Figure 3).

The following questions were further answered in this study. (1) How can orbital drivers (note that orbital obliquity was included while orbital eccentricity and precession were excluded) with limited periodicity produce period-abundant paleoclimatic variations within the PSMGIM? (2) Can the PSMGIM be significant for Earth’s heat budget? (3) How will PSMGIM variations contribute to paleoclimatic variations? The last question was partially explained and needed further studies.

2. Methods

This section contains three major parts, as follows: (1) the dynamic model for SMG-induced magma motion, (2) the periodicity calculations included for SMG-induced magma motion and the probability for a period to occur in SMG-induced magma motion, and (3) a spectrum analysis for data and modeling results.

2.1. The Dynamic Model for SMG-Induced Magma Motion

As far as methodology is concerned, omitting small-magnitude SMG in climate is an extension of the scale-analysis method that has been effectively applied for linear processes and short-term weather systems. However, for climate and paleoclimate studies with a large space and a long time duration, small driving factors such as SMG that accumulatively and actively acts on climate and paleoclimate systems may not be omitted. The size of SMG is even comparable to Coriolis in the “SMG dynamical zone (SMDZ)” [17], especially for slow magma. As compared to magma whose mass is approximately 3,400 times that of the atmosphere plus the oceans, the atmosphere and the oceans are only a thin layer of fluids having much smaller thermal and dynamic inertia or memory [17] to produce weather-climate variations driven by both SMG and solar radiation. However, magma, the third geophysical fluid, has been omitted in climate-paleoclimate studies, although magma has a large mass through which the thermal contribution can be significant even if solar radiation hardly reaches magma. SMG can drive magma to move and the kinetic energy from motion can be transferred into thermal energy via the friction associated with sticky magma.

Obtaining an accurate nonlinear solution for magma motion under SMG in Eulerian system is currently impossible. SMG changes with a relative location between a float and the Sun or the Moon. For a numerical model to determine the changing location of a moving float, grid spacing in Eulerian system must be smaller than the distance the float moves within one time step. The longer the temporal scale is required to study or predict, the smaller the speed dynamically contributes to the corresponding temporal variation [17], indicating that a higher resolution is required for a longer-term climatic model. Pinpointing accurate relative locations between a float and the Sun or the Moon given Earth’s, the fast Earth’s rotation makes the time step much shorter than that used by classic climate models. If, for example, the time step is one hour and the smallest speed that needs to be simulated is 10−5 m/s (as is typically required for a variation of ~1,000 years, Table 1), grid spacing must be smaller than 0.036 m. Increasing the time step may enlarge the grid spacing but will result in error for the relative locations and the momentum accumulation effect will be missed. Another limitation is due to the discretization requirement required for a stable and accurate numerical solution of determinacy (i.e., the well-posedness of a numerical solution), with a higher spatial resolution usually corresponding to a shorter time step [18].

To avoid difficulties associated with the Eulerian approach and to try to shed some light on SMG dynamic effects for paleoclimate, the dynamic equations of the SMG were established within the Lagrangian system via local coordinates on a rotating Earth. With a small Coriolis force, the widely used geostrophic balance approximation [19, 20] was applied in order to provide a circulation background (even if the horizontal pressure gradient force is not balanced with the Coriolis force, their balance can be omitted as compared to the much larger frictional force). Dissipation took the form of “Rayleigh friction” that fits low-speed fluids such as magma. After a simplification [17], the locations of magma, the Sun, and the Moon were expressed with the local motion speeds of fluids, the rotation speed of Earth, the orbital obliquities of the Sun and Moon, the Earth’s revolution speed, time, and the initial latitudes and longitudes of fluids and the Sun and the Moon in local coordinates. The applied dynamic equations for magma motion are as follows:where for the Sun and for the Moon. and indicate the zonal and meridional SMGIM components, respectively. and indicate the real zonal and meridional speed components, respectively, and maintain balance between the pressure gradient force and the Coriolis force. and are the density and pressure of the fluid, respectively. and are the geographical longitude and latitude, respectively. , , , and , for subscripts , respectively. . is the dissipation coefficient (). rad and rad (+/− in the Northern/Southern Hemisphere) are the obliquity of the orbital plane of the Moon and Sun, respectively. , , , , , , , , , , , , , , and for longitude and latitude (divided by time ) for a moving float at time , respectively. ,580,000 m and 3,260,000 m are the mass-weighed mean radii of the lower mantle and the outer core, respectively. m and m are the distance between the centers of the Earth and the Moon and Sun, respectively. rad/s is the mean apparent angular velocity between the Sun and the Earth. rad/s and rad/s are the mean revolution angular velocity of the Moon and Sun, respectively. (m3 kg−1 s−2) and (m3 kg−1 s−2), where, m3 kg−1 s−2 is the universal gravitational constant. kg and kg are the mass of the Moon and Sun, respectively.

For magma with a lower speed and a large friction, , the speed of magma and the periods of the present time step can be obtained using the speed () and the location of the previous time, with at least a 4th-order accuracy. The speeds (m/s), averaged from to seconds (or to days), were derived as follows:

Here, , , and .

The PSMGIM for horizontal motion can be written as follows:where kg and kg are the approximate mass of magma within the lower mantle and the outer core, respectively. The PSMGIM was directly computed from (1)–(3) using parameters located on the right-hand side of these equations.

2.2. Probability Included in SMG-Induced Magma Motion

The basic (named 1st, 2nd, and 3rd, resp.) period components included in (2) were the following:The comprehensive memorial periods can be quantitatively analyzed as follows:where .

For a period range , with a period , a mean length , and a corresponding speed range of , , the occurring probability (OP) for a specific PR within a total speed space (, , and , the horizontal components of the total speed) is designated as and derived in the and or directions, respectively, as follows:where and are the mean length of a -belt, the area circled by lines with a period of and within the speed space of , . and , the very small fractions, represent the area-overlapping percentage among various -belts determined using (5).

2.3. Spectrum Analysis of the Data and the Modeling Results

Spectrum methods [8, 9, 21] were used to determine the periodicity included in a time series. A squared amplitude spectrum method was developed in this study and used for a time series with a large sample size. The time series (here, temperature, insolation, CO2, dust, and the simulated SMG-induced power of magma) was at a discrete time , to Nt (where is the time interval of Pt and Nt is the sample size of Pt, yr, and ,000). The reconstruction of Pt via the Discrete Fourier Transform (DFT) is given by the following:

The discrete periods and the corresponding square amplitudes were, respectively, as follows:

Theoretically, the integer , but . Here .

(Some publications have not applied DFT correctly. As a result, the amplitude of time series reconstruction is much larger than that of the time series. Here, using a sample size of 8,000, time series reconstructions were almost identical to the time series and provided accurate spectra with a correlation coefficient of 100% and a confidence of 99.99% for all of the presented figures).

3. Results

3.1. The Significance of the PSMGIM

Magma is located within the lower mantle (at ~900 to ~2,900 km in depth and has a mass of ~2.94 × 1024 kg and a mass-weighed-mean radii of ~4,580 km) and the outer core (at ~2,900 to ~5,100 km in depth, with a mass of ~1.84 × 1024 kg and a mass-weighed-mean radii of ~3,260 km). Magma became the major fluid left within Earth as Earth cooled and collapsed with time and contributes to Earth’s heat budget via formation and disintegration [22–28] (e.g., Earth’s collapse, Earth’s inner core formation and separation, Earth’s mantle and core cooling, and Earth’s radioactivity within the upper mantle and crust). To maintain an energetic Earth, a reproducible energy source that does not originate from Earth’s formation and disintegration must exist and is likely not derived from solar radiation that hardly heats the inner Earth. The power and its variation for SMG-induced magma motion might help provide reproducible energy and explain paleoclimatic variations. The large contrast between magma and the atmosphere and oceans by mass (the mass of magma is ~80% of Earth’s mass, ~901,887 times the mass of the atmosphere, and ~3,400 times the mass of the oceans) makes PSMGIM a sizeable source for Earth’s heat budget, under gravitation from even larger celestial objects (here only the Sun and Moon were considered). Depending on the magma speed associated with dissipation (for given orbital properties and Earth’s constant rotation and radius), mean total PSMGIM could be ~800 to ~11,600 TW, with a standard deviation of ~500 to ~7,000 TW.

The “buffer” and “quick release” characteristics determined for the PSMGIM (with a smaller amplitude and a longer phase for negative anomaly phases as compared to positive anomaly phases) were derived from the experimental dissipation coefficient (no data available) set from 10−4 to 10−3 s−1 for magma within the Outer Core and enlarged by 100 to 200% for magma within the lower mantle (Figure 3). Smaller dissipation associated with magma with a higher temperature or with a smaller density was determined to induce faster magma motion and to produce a larger PSMGIM for shorter periods (e.g., 1-2, 20–25, and 40–45 ky). Larger dissipation associated with magma with a lower temperature or with a larger density tended to induce slower magma motion and produced a smaller PSMGIM during longer periods (e.g., 40–45, 70–80, and 100 ky). Additional experiments yielded a considerable PSMGIM. Approximately 800 to 11,600 TW for mean total PSMGIM (or ~500 to ~7,000 TW for its variation) is equivalent to ~1.57 to ~22.7 W (or ~0.98 to ~13.7 W) per square meter of Earth’s surface. The range is considerable when compared to the following: (1) the ~32 TW of heat flowing through the seafloor associated with thermal thinning of the lower lithosphere [29], (2) the W of power per square meter of Earth’s surface due to increased anthropogenic greenhouse gases and aerosols [30], and (3) the total heat budget of Earth (~41 TW) which includes Earth’s heat loss where approximately 54% cannot be balanced [22]. The PSMGIM may be one of the most reproducible external energy resources required to maintain an energetic Earth and to compensate for Earth’s heat loss.

3.2. Periodicity in the PSMGIM

SMG drives magma to move nonlinearly because the location of magma relative to the Sun or Moon changes with time. Momentum can be accumulated in moving magma under a rotating Earth on multiple periods that may be much more abundant than the original periods of the SMG, depending on orbital properties, including the following: (1) the revolution angular velocity of the Sun or Moon, (2) the relative velocity and the dissipation of magma motion, and (3) Earth’s radius and rotational velocity (as defined in (4)–(6) and simply demonstrated in the basic components of memorial periods in (4), as partially listed in Table 1). The occurring probability (OP) for a specific period range (PR) was quantified using (5)-(6) and is provided in Figure 2. In general, as shown in Figure 2, SMG-induced magma motion can have a broad-spectrum of periods that vary from 0 to infinity, with faster speeds associated with shorter periods, and vice versa, within the major period-speed area (beyond this area faster speeds may also lead to longer periods but with a much lower OP). OP was proportional to PR and Earth’s rotational speed (without rotation, OP would be zero for any period), but approximately inversely proportional to the total speed space and the square of the mean period length. Also, as indicated in (6), OP was larger at lower latitudes than at higher latitudes. For oscillation-preferring situations (e.g., at lower latitudes), the total speed space should enlarge. However, the enlarged total speed space decreased OP, indicating an “astringency” mechanism as for most natural systems. Figure 2 provides (6) using period ranges corresponding to the speed ranges for SMG-induced magma motion. Averaged from different latitudes, the OP of the PR was approximately as follows: (1) a PR of 1 to 3 ky had an OP of 59 to 69% for a speed space of −750 to 750 μm/s, (2) a PR of 15 to 30 ky had an OP of 48 to 56% for a speed space of −48 to 48 μm/s, (3) a PR of 35 to 60 ky had an OP of 44 to 52% for a speed space of −18 to 18 μm/s, and (4) a PR of 80 to 150 ky had an OP of 46 to 53% for a speed space of −9 to 9 μm/s.

Figure 2: Period ranges (PR, -labels) and their occurring-probability (OP, %). , set as the abscissa/ordinate, is the zonal/meridional speed component for magma within the lower mantle/outer core (μm/s, above/below the frame in Rows a, b, c, and d). Columns 1–3/4–6 provide the periodicity of zonal/meridional motion at the labeled latitudes (north is positive). OP is the ratio of the shaded area for PR to the total area labeled on the “” column, computed using (6). The OP averaged for Columns 1–6 is labeled on the “OP” column.

Figure 3: Squared-amplitude spectra for the total PSMGIM time series (TW2). The experimental dissipation coefficient for magma within the outer core was set at 10−4 (left column) and 10−3 s−1 (right column) and enlarged by 100% (Row 1), 150% (Row 2), and 200% (Row 3) for magma within the lower mantle. The total PSMGIM time series (not plotted) was 100% and was correlated with its reconstruction via the Discrete Fourier Transform (DFT) with a confidence of 99.99% and was evaluated using the following: (1) the maximum size of the positive/negative amplitude of its anomaly (“”), (2) the total phase length during positive-/negative-anomaly phases (“”), (3) the standard deviation (“STD”), and (4) the mean (“”). Parameters , , and are described in Section 2.

3.3. The Connection between the PSMGIM and Paleoclimatic Variations

What is the connection between the PSMGIM and paleoclimatic variations? Does the connection occurs through slow heat-transport or sporadic heat-eruptions associated with magma convection [26–31] from inside Earth? Or is it the result of the release of greenhouse gases (e.g., CO2) during periods of temperature changes for the solid Earth or during volcanoes [30] (~0.15 to ~0.26 gigatons of CO2 emissions via volcanoes per year)? Supposedly associated with eruptions, dust was determined to negatively correlate (−67.5%) with temperature and presented characteristics of longer accumulation and a quicker release during negative and positive anomaly phases, respectively, as compared to temperature (Figure 1(a) versus Figure 1(d)). Dust lagged temperature by approximately 500 years, implying that before eruptions that added dust to the atmosphere, heat from magma may have begun to heat the atmosphere via slow heat-transport and the release of greenhouse gases. Dust may have helped cool the atmosphere via negative feedbacks on temperature [32]. If only 1% of PSMGIM (e.g., 500 TW) were released, accumulating the energy (~1017 Joules) released during the eruption of Mount St. Helens in May 1980 [33] would only take ~6 hours. Therefore eruptions would occur too frequently due to the PSMGIM if the PSMGIM were not largely the result of slow heat-transport across a large area (~5.1 × 108 km2 at Earth’s surface). Such assumptions, however, are not yet proven.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.