Evidence of cosmic recurrent and lagged millennia-scale patterns and consequent forecasts : multi-scale responses of solar activity ( SA ) to planetary gravitational forcing ( PGF )

Solar activity (SA) oscillations over the past millennia are analyzed and extrapolated based on reconstructed solar-related records. Here, simple recurrent models of SA signal are applied and tested. The consequent results strongly suggest the following: (a) the existence of multi-millennial (∼ 9500-year) scale solar patterns linked with planetary gravitational forcing (PGF), and (b) their persistence, over at least the last glacial– interglacial cycle, but possibly since the Miocene (10.5 Myr ago). This empirical modeling of solar recurrent patterns has also provided a consequent multi-millennial-scale experimental forecast, suggesting a solar decreasing trend toward grand (super) minimum conditions for the upcoming period, AD 2050–2250 (AD 3750–4450). Taking into account the importance of these estimated SA scenarios, a comparison is made with other SA forecasts. In Appendixes A and B, we provide further verification, testing and analysis of solar recurrent patterns since geological eras, and their potential gravitational forcing.


Introduction
Long-term solar variability: It is possible that we have missed the forest for the trees.Driven by pragmatic hopes of finding keys to weather prediction, we run the risk of concentrating too much on timescales of practical consequences -of days, months or years.In taking a longer view we see the problem in clearer perspective. . .(Climate and the Role of the Sun; Eddy, 1981) Solar activity (SA) has non-linear characteristics that influence multiple scales in solar processes (Vlahos and Georgoulis, 2004).For instance, millennia-scale solar oscillations have been recently detected, like those of about 6000 and 2400 years, by Xapsos and Burke (2009) and Charvátová (2000), respectively, with important and interesting influences in the near, past and future climate.These millennial-scale patterns of reconstructed SA variability could justify epochs of low activity, such as the Maunder minimum, as well as epochs of enhanced activity, such as the current Modern Maximum, and the Medieval maximum in the 12th century.
Although the reason for these SA oscillations is unclear, it has been proposed that they are due to chaotic behavior of non-linear dynamo equations (Ruzmaikin, 1983), or stochastic instabilities forcing the solar dynamo, leading to on-off intermittency (Schmitt et al., 1996), or planetary gravitational forcing with recurrent multi-decadal, multi-centennial and longer patterns (Fairbridge and Sanders, 1987;Fairbridge and Shirley, 1987;Charvátová, 2000;Duhau and de Jager, 2010;Perry and Hsu, 2000).It should be noted that all proponents of planetary forcing have forecasted a solar grand minimum for the upcoming decades, but one of them has also forecasted a super minimum for the next centuries (Perry and Hsu, 2000).In addition, during recent decades, statisti-Published by Copernicus Publications on behalf of the European Geosciences Union.
cal forecasts (with physically based spectral information of reconstructed records) of solar magnetic activity predict a clear decrease in SA, reaching a minimum around AD 2100 (Steinhilber and Beer, 2013;hereafter S13;Velasco Herrera et al., 2015).
It should also be noted that several recent studies have been devoted to reconstruct multi-millennial-scale solar/climate related records.In relation to solar records, Steinhilber et al. (2009, hereafter S09;2012;hereafter S12), Solanki et al. (2004;hereafter S04), and Finkel and Nishiizumi (1997;hereafter FN97) have investigated isotopic concentrations in ice cores and tree rings over the past 9500, 9500, 11 500, and 40 000 years, respectively, in order to estimate SA and/or related isotopic production.Another two isotopic reconstructions by Sturevik-Storm et al. (2014;hereafter SS14) and Adolphi et al. (2014;hereafter A14) have just provided detailed key information over 20 kyr for the past interglacial, or Eemian, more than 100 000 years ago, and for the last deglaciation, over 8 kyr from 19 to 11 kyr ago, respectively.
These different cosmogenic radionuclide-based reconstructions of SA present variations for the past millennia, and as Muscheler and Heikkilä (2011) have pointed out, large uncertainties appear in reconstructions of the solar modulation of galactic cosmic rays from different proxies, 10 Be and 14 C, and of changes in the geomagnetic shielding influence.However, these reconstructed records provide, especially when considered all together, the most objective information as elements for detecting and eventually modeling and extrapolating multi-millennial-scale solar oscillations, trends and absolute levels.
In this work, we attempt to advance our knowledge of solar variability by considering reconstructed records of variables related with SA over the last glacial cycle, from isotopic information coming from ice cores and tree-ring layers, reanalyzing them with a linear modeling of oscillations with recurrent influences.This modeling is achieved through simple analogues and Fourier series models.Tests of the proposed method and the detected low-frequency solar signal, going back in time, are based on independent data.Finally, we discuss the different oscillations detected, the confidence of our forecasts, some alternative forecast methods, and astronomical information that suggests a possible planetary gravitational forcing of SA by an unknown mechanism during the last millennia.Further analysis is provided in Appendixes A and B. Appendix A: qualitative verification of total solar irradiance (TSI) ∼ 9.5 kyr recurrent patterns with a bivalve population (BVp) reconstructed from late Miocene (∼ 10.5 Myr ago) data; and Appendix B: empirical evidence of a lagged planetary gravitational forcing of the ∼ 9500-year total solar irradiance (TSI) recurrent pattern.

Data
In order to analyze SA recurrent oscillatory patterns, different reconstructed forcing records of these oscillations are used.
We have analyzed six different sets of solar-related information.Firstly, total solar irradiance (TSI or S, hereafter) reconstructed by S04, S09 and S12, based on the isotopic information of 14 C and 10 Be, have recently provided records of SA anomalies for the last millennia.Figure 1a displays S04, S09 and S12 reconstructed and intercalibrated values from 9450 BC to AD 1900, from 7360 BC to AD 2009, and from 7350 BC to AD 1988, respectively.The variance explanation (obtained as the square of the correlation coefficient multiplied by 100) between S04&S09, S04&S12, and S09&S12 for decadal average records are of 52.7, 82.9 and 59.3 %, respectively.
Secondly, there are three interesting and useful solarrelated records, that is 10 Be isotope concentrations, from Greenland ice cores: one covering the past 40 kyr (FN97), another covering only 20 kyr but belonging to the Eemian (SS14), and another covering 20-10 kyr BP (A14).Figure 1b and c display the information of 10 Be reconstructed records by FN97, SS14, and A14.

Modeling
To take into account different timescale recurrences, the solar/climate, SC, variable can be expressed with two models.One model is based on the Fourier Series (FS), another is based on a linear transformation of the proxy variable values, and the last one is based on temporal analogues.
The FS model can be written by means of Here, T is the FS base period, N FS represents the number of FS terms or harmonics, j is an index component term, a and b are amplitudes, t is time, and e FS (t) is the error in this model.
The analogue model is defined as Here, α A is the amplification factor, β A is the slope, δ A is the lag, γ A is the additive constant, t 1 is the initial time for the modeled period, and e A (t) is the analogue error of this model.In all these models, parameters are estimated through iterative or multi-linear regression processes that minimize the RMS values of errors.
Taking into account both the dating limitations and the approximated values provided by proxy reconstructions, and  (Finkel and Nishiizumi, 1997, FN97;Sturevik-Storm et al., 2014, SS14).(c) 10 Be isotope concentration in Greenland ice cores during the period from 19 to 11 kyr BP (Adolphi et al., 2014, A14).Please note that in all figures: as the 10 Be concentration varies inversely with solar activity, TSI, the beryllium scales are inverted, and thus the upper parts in these scales indicate high TSI levels.
instead of developing statistical analysis as convergence and confidence level estimations, we prefer in this stage of research on solar climate recurrences, to apply verification/replication of all of our findings with independent information in our estimation processes and results.Future climate reconstructions with more accurate information will provide further and refined statistical analysis.

Long-term solar-activity recurrent patterns
In order to detect multi-millennial-scale recurrences and/or persistent oscillations in SA, we need to analyze 10 Be information since it is a solar proxy variable and it is available over longer periods than SA records (SS14).However, there are several 10 Be post-production and fallout processes (i.e.residence time in the atmosphere, scavenging rate, troposphere-stratosphere exchange, precipitation rate, etc.) that may alter the concentration found in the ice archive (FN97, SS14).Accepting that 10 Be concentration variability is influenced by climatic variability through long-term variable trends and modulations, we propose to apply a homogenization process based on statistics to the 10 Be (FN97) record.Firstly, a detrending process based on polynomial expressions was applied.And secondly, a demodulation was applied in an attempt to make the variance uniform.The consequent results show the 10 Be atmospheric signal of this process with approximated recurrent oscillations with lags of 9.6 and 19.2 kyr, which are shown in Sect.S1 in the Supplement.
The statistically detrended 10 Be FN97 record was modeled with a periodic FS function with N FS = 10 that employed Eq. ( 1).After minimization processes, a 9390-year period, P , was found and the corresponding model that explain 49.2 % of variance is displayed in Fig. 2a.
It should be noted that the solar and climate recurrence periods evaluated with FS and analogue techniques (shown in Supplement) have shown values of 9500 ± 100 years (∼ 9.5 kyr).
Before extrapolating the 10 Be ∼ 9.5 kyr recurrences to TSI, we applied a wavelet analysis to the three TSI records.The TSI spectral results (see Sect.S4) show three main, significant periodicities around 5000, 2400 and 900 years, and confirm the existence of at least three harmonics of the ∼ 9.5 kyr oscillations in solar activity.
3.2 Verification of the recurrences of 10 Be ∼ 9.5 kyr patterns Although this FS periodic 10 Be model is based only on the last 40 kyr (see Fig. 2b), it was extrapolated to cover the last 130 kyr, for comparison with other independent information of 10 Be.A detailed comparison with the 10 Be SS14 record (in 5 parts) coming from Greenland and the Eemian is displayed in Fig. 2c.The maximum variance explanation, of 18.4 %, corresponds to a temporal adjustment of 2.5 kyr (a temporal bias going back in time) of the SS14 dating.This temporal adjustment is justified because a similar one, of 2.3 kyr, is required by the SS14 18 O record when it is compared with another reconstruction from NGRIP Greenland ice cores by Kindler et al. (2014;hereafter K14), which is shown in Sect.S2.This comparison constitutes an important verification and test of the proposed FS model.
In order to verify the detected recurrent patterns of 10 Be, we apply different homogenization and extrapolation processes to FN97 data.Specifically, we follow the original calculations made by FN97 and the suggestions provided by Dr. K. Nishiizumi (personal communication, 2014), and we have also calculated the atmospheric signal of 10 Be ( 10 Be Atm) based on accumulated snow (Cuffey and Clow, 1997) and the signal of 10 Be coming from the GISP2 ice core.Our normalizations, which are devoted to eliminating high-frequency local climate influences on the 10 Be signal, have provided elements (records) to confirm the previous results for the ∼ 9.5 kyr recurrence and a consequent increase and diminishing of the 10 Be Atm and TSI signals, respectively, for the following centuries, as also shown in the Sect.S3.
We have also shown, in Appendix A, a qualitative verification of the total solar irradiance (TSI) ∼ 9.5 kyr recurrent patterns with a bivalve population (BVp) reconstructed from late Miocene (∼ 10.5 Mya) data by Harzhauser et al. (2013).They developed a comparison, in the frequency domain, of their reconstruction, BVp, with solar activity, or TSI, (BVp, TSI) and found significant oscillations with periods from 20 to 200 years.Motivated by their findings, in an independent effort, we have developed the same comparison (BVp, TSI), but realized in the time domain over the complete 8000-year record, with an extrapolated and adjusted TSI pattern.Thus, our comparison was in a different domain, the time domain, which considered all the range of oscillations.

Empirical tests of the recurrence and potential
mechanisms of TSI ∼ 9.5 kyr patterns Looking for physical basis and robust evidence of the detected recurrences, we have developed two qualitatively different tests of the multi-millennial recurrences of TSI: a test based on a suggested gravitational forcing, and a test based on the TSI (S04) and 10 Be (A14) records.
The first test, based on a physical mechanism, which develops a gravitational forcing analysis, is shown in Ap- pendix B. It is an empirical analysis of gravitational forcing due to lateral forces.Those lateral forces generate a lowfrequency signal with a period of ∼ 9500 years, preceding by ∼ 6700 years, and is similar to low-frequency solar activity.Additional analysis of a non-linear lagged response of TSI to gravitational forcing is developed and suggests a logarithmic model variation for different gravitational forcing periods.
The second test, which is based on a high-resolution independent and normalized 10 Be record (A14), consists of an extrapolation backward in time of the TSI(S04) record (this record reconstructed sun-spot numbers or SSN), which is based on 14 C records.Although 14 C records are coming from well-dated tree-ring studies, a temporal bias correction of a 70-year lead was applied to the 14 C based TSI(S04) original record, before being extrapolated backward in time.This lead adjustment of 14 C records is justified, to compensate their limitations, because these records are influenced by the global carbon cycle, causing fluctuations of the atmospheric 14 C concentration measured as 14 C in tree rings to be damped, smoothed and delayed relative to the 14 C production (S12).
After this important adjustment, an application of the analogue model (a linear leaded transformation with corrected trend), Eq. ( 2), produces an excellent agreement between the 14 C based TSI(S04) record of SSN, with a lead of 9400 years, with TSI( 10 Be[A14]) record, displayed in Fig. 3.

Application of the ∼ 9500-year recurrence of SA
We applied Eq. ( 2) with a lag parameter of 9600 years to the TSI records, maximizing the match between the analogue model based on S04 information and the original S04 records.Only the S04 model continually covers the next centuries, due to its longest characteristics, and presents an overlapping that explains 16 and 53.4 % of the TSI variance of the last 1000 and 500 years, respectively.Results of TSI are displayed in Fig. 4. In this figure, the three TSI records (S04, S09 and S12) are displayed with their analogue models.Our model confirms a grand minimum in the period from AD 2050 to 2200 forecasted by S13, showing a sustained deficit of 0.5 W m −2 , similar to that shown in the Maunder minimum, 4 centuries ago (see Fig. 4b).
The same model, shown in Fig. 4a, suggests that the next super minimum of SA will occur from −2100 to −2600 BP (AD 4050-4550), and will be similar to the period 7500-7000 BP of reduced SA.Please note that in Fig. S1 in the Supplement, based only on 10 Be, an experimental forecasted decrease of TSI is shown with an SA super minimum from −2 to −2.5 kyr BP (AD 3950-4450).In Fig. 4, big and small vertical orange arrows indicate, the verified forecasts of super and grand solar minima, respectively.

Discussion
Firstly, and in order to confirm this multi-millennial recurrence, we have developed different tests and verifications of the SA recurrent patterns.In the following a summary of the tests and verifications of our findings is presented.
Our FS model explains the detrended and modulated 10 Be statistically corrected variability over almost the last 40 kyr.However the recurrent patterns based on FN97 when extrapolated backward in time are comparable with independent 10 Be information from the Eemian.
When this recurrent phenomena detected in the 10 Be record was extrapolated to the TSI records, we conducted other tests, establishing the following: (a) the overlapping of the TSI (S04) record explains over 53 % of the variance in the last 5 centuries; and (b) the extrapolated model also based on TSI (S04) presents an important match with different data (S12) and an independent procedure (FFT) employed in the TSI forecast due to S13.
In Appendix B, based on both the Solar System movement reconstruction and simulation, H services, developed by the Solar System Dynamics Group of the JPL/NASA, and monthly SSN data from the World Data Center SILSO (SILSO, 2015) for the last decades, we have provided additional elements to support the idea that long-term solar activity is modulated by recurrent planetary effects.Our analysis establishes the following: a. Solar System dynamics generate lateral forces (enhanced by its double integral) with multi-millennial- e.The similarity of the ∼ 9500-year TSI with the average SSN 10.5-year cycle, with scales differing at almost 3 orders of magnitude, suggests a self-similar process with a mechanism possibly linked to recurrent PGF in different scales.
The qualitative verification of the total solar irradiance (TSI) ∼ 9.5 kyr recurrent patterns with a bivalve population (BVp) reconstructed from late Miocene (∼ 10.5 Mya) data, shown in Appendix A, confirms not only the existence of this recurrent solar pattern throughout geological eras but also its persistent characteristics, period and pattern, because the comparison is made with the extrapolated forward in time TSI record (see Appendix A and Fig. 4).
Our experimental multi-millennial-scale analogue forecast of TSI, supported mainly by recurrent oscillations over the last glacial-interglacial cycle, shows a lowering trend toward a minimum for the coming decades.Our forecast also confirms previous efforts by several authors (Fairbridge and Sanders, 1987;Fairbridge and Shirley, 1987;Perry and Hsu, 2000;Duhau and de Jager, 2010), who have forecasted a solar grand minimum for the upcoming decades.For instance, recent findings linked to periodicities of the solar tachocline and their physical interpretation may permit us to estimate that solar variability is presently entering into a long grand minimum, thus consisting of an episode of very low SA (Duhau and de Jager, 2010).
Although the complete physical basis of this recurrent process is missing, there are several examples of physical and theoretical evidence that also support our findings.Firstly, it is important to highlight what Mackey (2007) has stated: "In several papers, Rhodes Fairbridge and co-authors described how the turning power of planets is strengthened or weakened by resonant effects between the planets, the sun and the sun's rotation about its axis."Specifically, there are important works motivated by Rhodes Fairbridge and other researchers, providing a theoretical basis and practical evidence of resonant interactions, for instance: a. Abreu et al. (2012) have shown the physical basis of a gravitational forcing of the solar tachocline variations.
They developed a gravitational model for describing the time-dependent torque exerted by the planets on a non-spherical tachocline and compared the corresponding power spectrum with that of the reconstructed SA record.They find an excellent agreement between the long-term cycles in proxies of SA and the periodicities in the planetary torque (with a period from 50 to 504 years).
b. Fairbridge and Sanders (1987) have indicated long-term variations due to planetary forcings.They follow Stacey (1963) who, based on the periodicities of planetary orbits, proposed a ∼ 4.45 kyr Outer Planets Restart (OPR) cycle.It is close to half of the ∼ 9.5 kyr detected periodic recurrence.
c. Looking for solar-planetary resonances of our detected ∼ 9.5 kyr, we compared the "biggest" solar system secular frequencies determined by Laskar et al. (2011) over 20 Ma for the four inner planets, and over 50 Ma for the five outer planets, corresponding to 45.184 and 49.880 kyr, respectively.We found that the mean value of 47.532 kyr is almost five times the solar period detected (47.532 kyr = 5× [9.56 kyr]).This means that the solar inner and outer planets show a resonance (5 : 1) with the solar periodicity detected.
We have found and tested a recurrence of ∼ 9500 years of SA that implies a solar grand minimum for the next one and a half centuries.However, we can also support our findings with other studies.For instance, the existence of different solar modes of activity (grand minima, regular, and a possible grand maxima), which have also shown important temporal variations with asymmetries (grand maxima significantly less often experienced than grand minima) during the Holocene (Usoskin et al., 2014), would be considered expressions of our detected recurrent pattern of ∼ 9500 years.
In this work, we have forecasted a continuation of the solar decline for the next decades, which is supported through precursory signals during recent decades: a.A steady and systematic decline in solar polar magnetic fields, starting from around 1995, which is well correlated with changes in meridional-flow speeds (Janardhan et al  (Janardhan et al., 2011).
c.A significant reduced ionospheric cut-off frequency to radio waves, normally about 30 MHz, to well below 10 MHz (Janardhan et al., 2015a).
Also, in this work, we have forecasted a grand solar minimum, with sustained low solar activity for the next 2 centuries, which has been supported through a number of recent studies and their findings: a.The continuation of this decline in solar activity is estimated to continue until at least 2020, and there is a good possibility of the onset of a grand solar minimum from solar-cycle 26 onwards (2031) (Janardhan et al., 2015b).
b. Based on the S04 SA record, it has been shown that gradual (abrupt) changes in solar surface meridional flow velocity lead to a gradual (abrupt) onset of grand minima, and that one or two solar cycles before the onset of grand minima, the cycle period tends to become longer (Choudhuri and Karak, 2012;Karak and Choudhuri, 2012).It is noteworthy that surface meridional flows over Cycle 23 (Hathaway and Rightmire, 2010) have shown gradual variations, and Cycle 24 started 1.3 years later than expected.

5
An analysis and test of recurrent solar variability for the last millennia has been presented in this study.It was based on multi-millennial solar-related reconstructed records from different and valuable proxy information.
The tested existence of the ∼ 9.5 kyr period recurrent pattern suggests that SA is characterized by solar dynamics with long-term patterns.Considering that it has been suggested that the modulating oscillations of SA, around 84, 178 and 2400 years, are possibly related to the Sun's rotation rate and impulses of the torque in the Sun's irregular motion (Landscheidt, 1999;Fairbridge and Sanders, 1987;Charvátová, 1995Charvátová, , 2000)), our results also suggest that similar mechanisms on the solar dynamo must be proposed for solar oscillations of around 9.5 kyr.This hypothesis should be tested, taking into account the results presented in this paper.
In this direction, we analyze two key evidences of the solar recurrent pattern.In Appendix A, we present a qualitative verification of the total solar irradiance (TSI) ∼ 9.5 kyr recurrent patterns with a bivalve population (BVp) reconstructed from late Miocene (∼ 10.5 Mya) data that show the geological-scale persistence and regularities of the SA patterns, which can only be explained by a planetary gravitational forcing (PGF).
In Appendix B, we present an empirical analysis of solar PGF forcing due to lateral forces.We found that lateral forces generate a low-frequency (∼ 9500 years) signal that presents similarity with, and precedes by ∼ 6700 years, low-frequency solar activity.This lag appears to be part of the non-linear lagged responses of solar activity to different time-lengths of PGF.We suppose that this lateral forcing could enhance the oblateness of the solar body, and tidal influences, and consequently, as Abreu et al. (2012) have also suggested, regular cycles of solar activity.Finally, we also find that the solar patterns of 9600 and 10.5 years are similar, suggesting a common gravitational origin.
With all of these recurrent and persistent phenomena, we have presented, tested and verified an experimental multimillennial forecast technique for SA.We have provided elements and recent supporting studies on precursor signals of an entering into a grand minimum SA mode.The extreme duration of the last solar minimum is important evidence of longer cycles, similar to those presented before the start of the Maunder and Spörer minima.
We can conclude that the evidence provided is sufficient to justify a complete updating and reviewing of present climate models to better consider these detected natural recurrences and lags in solar processes.this forcing have been analyzed by only a few of them.For instance, Abreu et al. (2012) analyzed the PGF of solar tides.In this appendix, we compare the lateral (perpendicular to movement) forces evaluated by lateral accelerations of the solar movement around the solar system (SS) barycenter (BC) with TSI reconstructed and extrapolated (based on the detected recurrences) records.Our work is motivated by the findings from Fairbridge and Shirley (1987), who predicted the initiation of a Maunder-type prolonged solar minimum on the basis of a study of solar motion with respect to the SS-BC in the years from AD 760 to 2100.Their study detected "patterns" in solar orbits associated with different levels of SA.
In order to analyze variability in both solar dynamics and solar activity, we studied data coming from the following: (a) lateral forces (F ) of the sun due to planetary gravitational forces and movements, reconstructed and forecasted by JPL (2013)/NASA from 3000 BC to AD 3000, and (b) solar activity expressed in the total solar irradiance (TSI) average from three reconstructed records over the last millennia and extrapolated for the next millennia (all shown in Fig. 4a).These two records are displayed in Fig. B1.Lateral (perpendicular to movement) forces were evaluated based on x, y and z coordinates and derivatives provided by the HORIZONS (H) system from the Jet Propulsion Laboratory/NASA (JPL/NASA) for the past 5000 and future 1000 years, every 90 days, where the x-y plane is the ecliptic plane centered on the BC.As an approximation, solar movement was considered only in the x-y plane, and lateral acceleration in this plane was evaluated to be perpendicular to the tangential direction of solar movement.These on-line solar system data and ephemeris computation services provide accurate ephemerides for solar system objects.
The simulated lateral inertial forces (F ) are considered to provide gravitational influences on solar activity.In order to enhance their low-frequency oscillations, we applied the double integral function to the analyzed F record.
An integration was applied twice to the Solar signal, S(t), as follows: where S(t) is the solar signal, expressed as force (F (t)) or as total irradiance (T (t)), t is time, t 0 is initial time, σ S N (t) is its time integral, N is the successive application number, and µ S and µ S 1 are the long-term averages of S(t) and σ S 1 (t), respectively.Vertical scale of values were inverted both in σ F 2 (t) and σ T 2 (t) because the double integral procedure changes the sign of the enhanced result.Please note that the last minimum of lateral inertial force F was around 2500 BC and the next minimum of TSI will be expected around AD 4200.Please note the time difference between these minima of around 6700 years.
We apply Eqs.(B1) and (B2) to these two records in order to enhance low-frequency variations.Results are displayed in Fig. B2.The double integral of forces F , σ F 2 (t), is almost explained (99.9 % of variance) by a sine function of a 9400-year oscillation, which is also displayed in Fig. B2a.The double integrated solar activity TSI, σ T 2 (t), also shows a periodicity of ∼ 9500 years.The scales of these enhanced solar signals are inverted because the sign is changed due to the double integration enhancement.A comparison of both curves is displayed in Fig. B3a. Figure B3b also displays both integrated curves, however the σ F 2 (t) curve is lagged (moved forward in time) 6700 years.
In order to verify this 6700-year lag of the TSI response to the F oscillation of a 9500-year period, we look for two other similar pairs of periods and lags (P /L).In order to obtain an additional pair of P /L, we analyze the lateral force F and solar activity TSI over the AD 1000-3000 period.We applied a double integration (Eqs.B1 and B2) and a polynomial detrending process to F .The σ F 2 (t) and its trend are shown in Fig. B4.The detrended σ F 2 (t) is compared with the TSI record, and oscillations of ∼ 950 years are detected, together with a lag of ∼ 350 years of TSI with respect to the supposed forcing F .These two variables are displayed in Fig. B5.
Another pair of P /L values is evaluated with the Hale SSN cycle of ∼ 22 years that shows an alternating magnetic sign for each 11 years.It is compared with the Fourier series (with only two harmonics) of the force F , signal based on www.earth-syst-dynam.net/7/583/2016/Earth Syst.Dynam., 7, 583-595, 2016  the period from AD 1700 to 2000.The comparison, depicted in Fig. B6, indicates a lag of less than 1 year.The three sets of L/P pairs are 0.5/22, 350/950 and 6700/9500 years, respectively.These three pairs, which correspond to different phases, are modeled together with a nonlinear function that tends toward a lower limit of 0 • for lower periods, and an upper asymptotic limit of 360 • (2 * Pi radians).The adjusted model for phase variations in terms of period, which is a logarithmic function, is depicted in Fig. B7.
It is important to emphasize that the results in Fig. B5 have not only detected an L/P pair but also confirm the next forecasted grand minima (AD 2020-2220) also associated with contributions from ∼ 350-year lagged influences of solar lateral forces.Additionally, we developed an interesting comparison that shows self-similarity in TSI.This comparison is for our enhanced σ T 2 (t) ∼ 9500-year-solar-cycle, evaluated previously, with a Fourier series model of the solar SSN cycle with a period of 10.5 years, based on monthly SSN data from the World Data Center SILSO, Royal Observatory of Belgium, Brussels, over the period from 1964 to 2008.This comparison is displayed in Fig. B8, and clearly shows selfsimilarities between these two solar cycles.Both cycles show a shorter increasing period (∼ 25 %) than the decreasing period (∼ 45 %), and a maximum plateau (∼ 20 %), and an almost nonexistent minimum plateau.
Finally, we also developed a spectral analysis of the analyzed lateral forces.This analysis is based on wavelets and is displayed in Fig. B9.It clearly shows important contributions to periods around of 12, 22, 60 (in a range of 50-80), 180, 650 (in a range of 400-800), 1000 and 2600 years.

Figure 1 .
Figure 1.Solar-related signals.(a) Solar activity, TSI, reconstructed by S04, S09 and S12, after an intercalibration using the S09 record as a base.(b) 10 Be isotope concentration in polar ice cores during the past 130 kyr (Finkel and Nishiizumi, 1997, FN97; Sturevik-Storm et al., 2014, SS14).(c) 10 Be isotope concentration in Greenland ice cores during the period from 19 to 11 kyr BP(Adolphi et al., 2014, A14).Please note that in all figures: as the 10 Be concentration varies inversely with solar activity, TSI, the beryllium scales are inverted, and thus the upper parts in these scales indicate high TSI levels.

Figure 2 .
Figure 2. Data and modeling of 10 Be isotope concentration in Greenland ice cores.(a) Data and an FS model of 10 Be isotope for the past 40 kyr provided by Finkel and Nishiizumi (1997, FN97) after being detrended and demodulated (see Supplement).(b) The model, shown in (a), is extrapolated covering the last 135 kyr, and the (Sturevik-Storm et al., 2014, SS14) data are included for comparison.(c) A zoom of (b) for a detailed comparison of the extrapolated FS model and the SS14 data.Please note that a maximum match implies a SS14 temporal adjustment, or time bias, of 2.5 kyr going back in time.

Figure 3 .
Figure 3.A test of the recurrent TSI signal based on 14 C over the last 11 000 years, TSI(S04), and the linear model based on the 10 Be isotope concentration record from Greenland ice cores.We extrapolated the TSI(S04) record backward in time, 9400 years, to match the model based on 10 Be isotope (A14) that covers the period from 18 000 to 10 000 years BP.

Figure 4 .
Figure 4. Solar activity signals reconstructed and modeled records.(a) Solar activity, TSI, reconstructed by Solanki et al. (2004), Steinhilber et al. (2009, 2012) (S04, S09 and S12, respectively), shown in Fig. 1a, and their analogue models.(b) A zoom of (a) that covers only 2 kyr.(c) Another greater zoom of (a) that covers only 0.8 kyr including the independent TSI forecast by S13.Big and small vertical arrows indicate super and grand solar minima.

Figure A2 .
Figure A2.Comparison of bivalve population (BVp) and TSI adjusted with a linear transformation [TSIadj = TSI +a(t − t 1 ) + b] with a = 0.03 [W m −2 ]/1000[yr] and b = −0.3[W m −2 ].Only positive TSI adjusted anomalies are displayed in order to enhance the match and show a possible threshold provided by TSI for BVp development.

Figure B1 .
Figure B1.(a) Lateral forces on the Sun generated by planetary gravitational forces (PGF), expressed in [arbitrary units], evaluated by the Horizon/NASA system from 3000 BC to AD 3000.(b) The solar activity (TSI) expressed in [W m −2 ], average values of the three reconstructed records over the last millennia and extrapolated for the next millennia (10 000 BC to AD 12 000) based on ∼ 9500year recurrence (records shown in Fig. 4a).

Figure B2 .
Figure B2.(a) Double integral of the solar lateral inertial forces σ F 2 (t) due to the solar movement resulting from the planetary gravitational forces shown in Fig. 1a, and (b) double integral of the solar reconstructed and extrapolated record σ T 2 (t) shown in Fig. 1b.Vertical scale of values were inverted both in σ F 2 (t) and σ T 2 (t) because the double integral procedure changes the sign of the enhanced result.Please note that the last minimum of lateral inertial force F was around 2500 BC and the next minimum of TSI will be expected around AD 4200.Please note the time difference between these minima of around 6700 years.

Figure B3 .
Figure B3.A comparison of the double integral of the solar reconstructed and extrapolated record σ T 2 (t) and the model based on solar lateral inertial forces, σ F 2 (t), σ T 2 (t) [σ F 2 (t)].To enhance a possible cause-effect relationship [σ F 2 (t)-σ T 2 (t)] the σ T 2 (t) [σ F 2 (t)] record is shown (a) without and (b) with a lag of 6700 years.Vertical scale of values were inverted because the double integral procedure changes the sign of the enhanced result, thus upper/lower values indicate maxima/minima.

Figure B4 .
Figure B4.The double integral of the solar lateral inertial forces σ F 2 (t) for the period from AD 1000 to 3000.A polynomial trend is also depicted.

Figure B5 .
Figure B5.A comparison of the double integral of the solar lateral inertial forces σ F 2 (t) and the solar reconstructed and extrapolated record TSI.To enhance trends, polynomials are adjusted to both records, showing oscillations of ∼ 950 years and a lag of ∼ 350 years.

Figure B6 .
Figure B6.A comparison of the Hale solar cycle of SSN and the Fourier Series (with only two harmonics) model of the F signal.A lead of 0.5 years is applied to the SSN to improve the match with the F recurrent model.

Figure B7 .
Figure B7.A model of phase between F (forcing) and TSI (suggested response) signals for different periods.The adjusted logarithmic model, which explains 99.5 % of the variance, shows two trends, one to zero for short periods, and other, an asymptotic one, for long periods.

Figure B8 .
Figure B8.Comparison of (a) the double integral of the reconstructed solar activity (TSI) (Vertical scale of values were inverted because the double integral procedure changes the sign of the enhanced result), and (b) the mean solar cycle obtained with an FS model based on SSN data from the World Data Center SILSO, Royal Observatory of Belgium, Brussels, over the period from 1964 to 2008.

Figure B9 .
Figure B9.(a) Lateral force [arbitrary units].(b) The wavelet power spectrum.The contour levels are chosen so that 75, 50, 25, and 5 % of the wavelet power is above each level, respectively.The cross-hatched region is the cone of influence, where zero padding has reduced the variance.(c) The global wavelet power spectrum (black line).The dashed line is the significance for the global wavelet spectrum, assuming the same significance level and background spectrum as in (b).Reference: Torrence and Compo (1998).