A SIXTEEN-YEARS DELAY IN SUN-RELATED TERRESTRIAL WARMING

. In an earlier paper De Jager et al. (2018) studied the relation between the intensity variations of the variable equatorial and polar solar magnetic energy sources and that of the average terrestrial northern hemisphere ground temperature. It was found that during the period 1610 to about 1920 the variation in the smoothed terrestrial temperature and those of the solar magnetic radiation are well correlated but that after that an additional temperature increase became prominent, to reach around the year 2000 a value of about one degree above the sun-dependent forecast. In the present paper an additional aspect of the above material is described, and we find that, apart from the above mentioned features, there is a delayed component of the sun-related heating with a time lag of about 16 years. The resulting relation between the variable flux of the solar magnetic energy and the average - smoothed – northern hemisphere ground temperature is presented in our final diagram, Fig.13. Our results confirm earlier findings by other authors that ascribe such a time lag to glacier warming and cooling.


The solar influence on terrestrial temperatures.
We refer to the study by De Jager et al. (2018). In this investigation, based on the new sunspot counting system introduced by Clette et al (2012Clette et al ( , 2014, it was found that, during the period 1610 (first telescopic sunspot observations) till about 1920, the smoothed terrestrial temperature variations were well correlated with the solar variations of equatorial and polar activity. Both equatorial as well as polar variations affect the terrestrial temperature variations, but after these years the earth's northern hemisphere average ground temperature rose to reach a value of ~ one degree above the sun-related value in 2000.

Previously derived solar temperature models
A short introduction of the previously derived models in De Jager et al. (2018) is needed, because they are used as input for this study. During the period 1610 till 1920 the fair agreement between 2 solar heating and the earth's average northern hemisphere ground temperature showed a remaining residual curve with a standard deviation of 0.12 K. In the attempts to explain these additional residuals, the search was directed towards finding a significant correlation as a function of time, t, between the two variables: solar variability and terrestrial ground temperature variation. The resulting input data streams for the period 1610 -1910, taken from de Jager et al. (2018), are given in Fig. 1. Here, the dark curve shows the smoothed variation of the average NH ground temperature, that we call the MobBroh curve (cf. De Jager et al. 2013), while the green curve, this being the solar model GSN-1612, represents the smoothed values of the solar surface magnetic fields.  Table 1 shows all found solar temperature models in de Jager et al. (2018), which are based on the maximum of the group sunspot numbers, max Gn , and the minimum of the magnetic antipodal amplitude, min aa . The model names consist of the abbreviated proxy type with the starting year of the data set. An extra letter indicates a single proxy model for models with 2 proxies. The models in Table 1 are obtained by the least square method of curve fitting (Kutner et al., 2005; chapter 6 and 7). The models are described with a linear dependency on the proxies:

Solar Model Name Fitted coefficients Quality of fit
Columns 2 till 4 in Table 1 show the fitted coefficients a, b and const with their standard deviation between brackets. In columns 5 till 8 quality figures of the regressions are shown. Columns 5 till 7 show the t-test value of the fitted coefficients, which is the ratio of the mean value and the standard deviation of the fitted coefficient. An absolute t-value of 2.6 indicates 99% confidence that the value is different from zero, and a higher t-value gives more confidence. The last column shows the adjusted coefficient of determination 2 adj R , which indicates the fraction of the variance that is explained by the model. More detailed information about the used multiple regression method and explanation of the Table, can be found in de Jager et al. (2018). The models in Table 1 are used to determine the time lag at the highest correlation with the MobBroh ground temperature, and thereafter the improved solar temperature regression models are calculated with this optimal time lag.

Methods -Calculation of time lag in the terrestrial response with the GSN1612 model
The following analyses are based on the GSN-1612 solar model. The purpose is to correlate two data streams with data at one year time intervals over a limited time period with use of the Fourier pair of the correlation theorem (Press et al., 1992; chapter 12.0 and chapter 13.2). In the time domain the cross correlation curve of two data streams is found by shifting the streams past each other and summing the multiplications of each corresponding data stream point k with each other at each shift step j , as shown by Eq. (2): On the other hand correlation can be written as a convolution. Convolution is calculated by the same process as described for the correlation, but with one data stream reversed. Convolution or filtering of the two streams with each other is a simple multiplication process of each corresponding data point in the frequency domain. The reversed stream can be noted by a complex conjugation. The Fourier pair for cross correlation is: Before performing the multiplication in the frequency domain, both data streams will be transformed by an FFT (Fast Fourier Transform) and the result of the multiplication will be transformed back to the time domain by an IFFT (inverse FFT). Some data processing has to be done to get a good transformed result; viz. a 'detrended' one to get the variations 'above' the DC value and zero padding to interpolate in the frequency domain for more 4 accurate defined peaks. The detrend is made by a linear fit through the data stream and by subtracting this fit from the data stream. Fig. 2 shows the detrended data streams with slopes and offsets removed.

Figure 2 -Detrended input data streams; the slopes and offsets are removed
The data streams in Fig. 2 have 299 points. The FFT routine requires a n-power of twice the number of data points and so the amount of points must be increased. Hence, 2048 points were chosen by adding zeros (zero padding) around the data streams. The zero padding will give spectral interpolation by a factor of about seven in the frequency domain. Thereafter both streams are transformed by an FFT routine and the cross correlation spectrum is calculated by multiplying the frequency curves with each other (the MobBroh spectrum is complex conjugated). Fig. 3 shows the transformed streams and the resulting cross correlation spectrum.

Figure 3 -Both streams transformed to the frequency domain and the calculated correlation spectrum
The cross-correlation spectrum is then transformed back to the time domain with an IFFT function. If interpolation is desired in the time domain, zero padding can be performed at the short period sides of the spectra. With a double-sided spectrum, the same amount of points (again an n-power of 2) has 5 to be added to the negative frequency side. After performing a back transformation the cross correlation curve in the time domain is obtained. It was the aim to obtain the normalized cross-correlation or the correlation coefficient, and hence the time domain cross correlation curve should be scaled with regard to the maximum value. This maximum value can be found by the Cauchy-Schwarz inequality and is: The maximum expected value is the root out of the product of the auto-correlation of each data stream at zero lag. The auto-correlation at zero lag is the energy of the signal. The energies in the time domain are calculated from the data streams in Fig. 2 Fig. 4 shows the normalized cross correlation curve as a function of the lag. The correlation coefficient has a range of -1 to +1. Fig. 4 shows that the maximum cross correlation is not at zero years lag. and 30 years. It appears that at a lag of zero years the terrestrial response is immediate (with no data shifting being applied). Then, up to a delay period of 16 years it grows to a maximum, where the response of the earth has about 5% larger correlation magnitude compared to the value with zero lag, after which the magnitude diminishes.

The time lag in the terrestrial response with the GSN-1844 models
The same calculation steps are applied to the three GSN-1844 models. Fig. 6 shows the results with solar model GSN-1844. When considering the detrended data in the upper panel, it appears that the lag has disappeared. The reason can be found in the two proxy models: while both input streams of the model fit have about the same lag, what can be found in the single proxy Figs. 7c and 8c, the proxies have opposite signs. The delay information is cancelled by these opposite signs. While the fitted model has a better quality, the quality of the delay information is degraded and calculating the lag with this two proxy model is therefore not correct. A better model approach may be to calculate the lag of the proxies separately, and then refit the data with the two lag adjusted proxies; this method has not been applied yet.
As next steps the lag calculation is applied at solar model GSN-1844b with the maximum values of the sunspot group numbers max Gn , and at solar model GSN-1844c with the minimum values of the magnetic antipodal amplitude min aa . Note, that while GSN-1844c is a model with a single magnetic proxy, a better model name would be AA-1844 (where AA refers to the antipodal amplitude).
The results of the lag calculation with solar model GSN-1844b is shown in Fig. 7. The upper figure shows clearly that there is a lag present in the detrended data; the resulting lag is 7 years. When this result is compared with the result of the GSN-1612 model, it appears that it is not the same value. This leads to the assumption that the lag may not be constant in time. This conclusion is further examined in chapter 4.  Table 2 recaps the lag results for all models.

Solar Model Name Lag [years]
Corresponding graphs GSN

Results -Improved solar temperature models
The solar models can be improved by using the found lags of Table 2. This is done by shifting the proxy data further in time with the found lag, after which regression is performed with the MobBroh NH ground temperature data and the corresponding proxies. A detailed description of the applied multiple regression method can be found in de Jager et al. (2018). The improved model of GSN-1612, now called GSN-1628-S and in which the "S" indicates a model with proxy shifting, the max Gn data stream is shifted 16 years forward in time. This causes a loss of 16 degrees of freedom, because the measurements above 1910 cannot be used. The reason is, as was proven in De Jager et al. (2018), that the "modern temperature increase" starts after about 1915; these points will degrade the model.  Another attempt was made to fit GSN-1853-S with both shifted proxies, but it fails with a bad t-test (see Table 3). The reason is that both proxies are too well correlated; our multiple regression method cannot handle this correlation. A solution can be Ridge Regression or a linear Bayesian method with an appropriate prior distribution, but this has not been investigated yet. The resulting model parameters are shown in Table 3, which is the updated version of Table 1. It is obtained with least square regression at maximal correlation of the data streams. Bold number with an asterisk are failed reliability tests, which indicates a bad model solution.

Fitted coefficients Quality of fit
When Table 3 is compared with Table 1, it appears that the standard deviations are slightly better. The same conclusion can also be drawn from the higher t-test values: lower standard deviations with about the same mean value. The explained variance by the models has raised dramatically, which indicate more preferred solutions.

Preliminary discussion on variable lag between NH temperature and solar models
In the previous chapter it appeared that the GSN-1612 model lag is 16 years, but the GSN-1844b model resulted in a 7 years lag. A possible explanation is that the lag may not be constant in time.
This assumption is next examined with an automated script that calculates the lag according to the method used in chapter 2, but now in 35 years windows. The 35 years window is shifted up by 1 year after each lag calculation until the end of the dataset has been reached. Zero padding is applied in the frequency domain to get a more accurate value of the lag in the time domain. The result and the data trend, a second order exponential fit, are shown in Fig. 12.
The earlier found lag of the GSN-1612 model is roughly the average of the red points in Fig. 12 and the same applies to the GSN-1844b lag. We conclude that the lag is not constant in time. Scafetta et al. (2016) describe that the cycles in the 14 C  cosmogenic radioisotope record are related with the planetary mass center relative (PMC) to the sun. This PMC travels in loops around and through the sun, according to the positions of the planets. The lag variations in Fig. 12 are in the apo-/pericycle duration range of 7 till 16 years, also mentioned by Scafetta et al. However our data 13 streams are low pass filtered by the Lowess algorithm for periods shorter than about 18 years; hence only the Jovian planet resonances and subharmonics are present. Beer et al (2018) make the connection between the cycles in the 14 C  record and the cycles in the modulation of the Sun's magnetic field. Yndestad et al (2017) have studied the connection between the solar system oscillation (mutual gravity attraction between the sun and the planets) and the Total Solar Irradiance (TSI). In Yndestad et al. (2017) the frequency content of sun-related data sets, also including a sunspot data series, has been determined by wavelet analysis. This result is added as a sinewave with only the phase information of de Vries/Suess cycle, which is a 5/2 subharmonic related to Uranus with a period of 210 years, at the bottom of Fig. 12. It appears that the peak coincides with the peak of the lag trend curve. When also the phase of the 9/2 subharmonic from the TSI-LS data set from Yndestad et al. (2017) is added to Fig. 12, it can be seen that the phase of the lag trend cycle seems to be in anti-phase with the 9/2 subharmonic with roughly the same period of 378 year. This 9/2 subharmonic is also mentioned in Duhau et al. (2020) as 2-Suess or the double Suess mode with a duration in this particular period of 419.8 years. Of course one should be aware of the limited length in time of the present data set; more data is needed for a more solid conclusion. This topic will be continued in our next study.

Conclusions
The main results are summarized in Fig. 13.
We thus conclude that there is a significant time lag of some 16 years between solar heating and part of the terrestrial response. This exercise is not in disagreement with a result found earlier by Weiss (2010, see in particular his figure 13 and the relevant text). From a study of ice cores from the Belukha glacier in the Altai mountains of Central Asia a shift of about 20 years was found. The description of his figure 13 states that a correspondence between the ice temperature record and the solar modulation potential was found "after allowing for a lag of 20 years". This time lag is of the same order as the shift found here. These findings present results that at the same time confront us with some problems that, however, are mainly in the meteorological domain: -Can we find a mechanism that explains a time delay of some 20 to 30 years between solar input and terrestrial response? -If this time delay is real, can it be assumed to be responsible for (part of) the extreme warming around 1990 -2000 in relation to the exceptionally large solar activity peak that occurred near 1960? Some suggestions for a positive answer to these questions have been forwarded by Bucha (1983) and summarized by Duhau and De Jager (2016), from whom we quote: Evidences has been found by Bucha (1983) that when the geomagnetic field is mainly dipolar, and the geomagnetic poles are near to the geographical poles, as it has been most of the time during the last 50.000 yrs., solar storms may contribute to heat the atmosphere at the auroras and the polar cap, leading to changes in global circulation (Bucha, 1983). It also substantially contributed to melting the ice at higher latitudes, with increasing ice-albedo feedback of insolation (Duhau and Martinez 2012), leading to cloud cover variations ... that also may contribute to heat the higher latitudes by cloud cover feedback of solar irradiance, as model computations indicate (Holland & Bitz, 2003).
Hence there are indications that the time delay of 16 years, found in this study, and that of about 20 years, found by Weiss (2010), between solar input and terrestrial response, may be related to the time that it takes for relatively stronger melting of the major terrestrial ice masses. In agreement with this it was found that the glacier retreat after the Little lce Age is occurring since around the years 1800-1850 and has not ended yet (Masiokas et al. 2010;Akasofu, 2010). This melting of high latitude ice leads to a non-linear increase of the ice-albedo feedback of solar 15 irradiance, which in turn leads also to cloud cover feedback and thus may be contributing to the modern temperature increase (Again, this problem is mainly for the field of climatology and will not be dealt with here). -Total Solar Irradiance -error mean square -error sum of squares -regression sum of squares -total sum of squares -degrees of freedom -years

Availability of data and materials
The used data sets will be made accessible at www.cdejager.com/sun-earth-publications/ The Sunspot Group Numbers data sets are available at www.sidc.be/silso/datafiles