Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Solar-climate connections over the Common Era

Authors

Julien Emile-Geay ORCID, Feng Zhu ORCID and Deborah Khider ORCID

Preamble

Variations in solar irradiance have been hypothesized to be responsible for part of the variability in global mean surface temperature (GMST) over the Common Era (the last 2,000 years or so). In this notebook we use Pyleoclim to test this hypothesis using wavelet transform coherency between reconstructions of solar irradiance and GMST.

This notebook reproduces Fig S16 of PAGES2k Consortium (2019) and adds an analysis of phase statistics.

Technical Skills Involved:

Data

Reading time

5 min

Keyword

Pyleoclim; Coherence Analyis; Common Era; Solar irradiance

Let’s load the necessary packages:

# load packages
import pyleoclim as pyleo
import pandas as pd
import numpy as np
import os
from matplotlib import gridspec
import matplotlib.pyplot as plt
from tqdm import tqdm

Exploratory Data Analysis

Let’s start by opening the forcing data and GMST reconstruction:

database_path = '../data/'
forcing_filename = 'Neukom2019_forcing.csv'
GMST_filename = 'Neukom2019_GMST_recon_medians_selection.txt'

forcing_df = pd.read_csv(os.path.join(database_path, forcing_filename)) 
GMST_df = pd.read_table(os.path.join(database_path, GMST_filename))
display(GMST_df.head())
Loading...
display(forcing_df.head())
Loading...

Let’s import this information into Pyleoclim and plot the corresponding timeseries:

ts_solar = pyleo.Series(time=GMST_df['Year'].values, 
                        value=forcing_df['solar'].values, 
                        label='Total Solar Irradiance', 
                        time_name='Year', 
                        time_unit='CE', 
                        value_name='TSI', 
                        value_unit=r'W m$^{-2}$')

Time axis values sorted in ascending order

exp_names = list(GMST_df.columns.values)[1:]

ts = {}
for i, exp in enumerate(exp_names):
    ts[exp] = pyleo.Series(time=GMST_df['Year'].values, 
                           value=GMST_df[exp].values, 
                           label=exp, 
                           time_name='Year', 
                           time_unit='CE', 
                           value_name='T anomaly', 
                           value_unit='K')
Time axis values sorted in ascending order
Time axis values sorted in ascending order
Time axis values sorted in ascending order
Time axis values sorted in ascending order
Time axis values sorted in ascending order
Time axis values sorted in ascending order
Time axis values sorted in ascending order

pyleo.set_style('web')
gs = gridspec.GridSpec(4, 1)
gs.update(wspace=0, hspace=0)

fig = plt.figure(figsize=[10, 6])
ax1 = plt.subplot(gs[0:3, :])
for i, exp in enumerate(exp_names):
    ts[exp].plot(ax=ax1,linewidth=1)
ax1.legend(bbox_to_anchor=(1.3, 1.0), loc='upper right')
    
ax2 = plt.subplot(gs[3, :], sharex=ax1)
ts_solar.plot(ax=ax2, linewidth=2, color='black')    
ax2.legend(bbox_to_anchor=(1.21, 1.0), loc='upper right')
<Figure size 1000x600 with 2 Axes>

Wavelet coherency analysis

This uses Pyleoclim’s wavelet coherence method. Since the series are evenly-spaced, we will use the CWT method.

scale = np.logspace(1, 11, num=51, base=2)
freq = 1/scale[::-1]
coh = {}
for exp in tqdm(exp_names):
    print('Evaluating coherence with '+ exp +'\n')
    coh[exp] = ts_solar.wavelet_coherence(ts[exp], 
                settings={'freq':freq}).signif_test(number=200, mute_pbar=True)
  0%|          | 0/7 [00:00<?, ?it/s]
Evaluating coherence with CPS_R-FDR_longcalib


 14%|█▍        | 1/7 [00:06<00:41,  6.95s/it]
Evaluating coherence with PCR_R-FDR_longcalib


 29%|██▊       | 2/7 [00:14<00:35,  7.04s/it]
Evaluating coherence with M08_R-FDR_longcalib


 43%|████▎     | 3/7 [00:21<00:28,  7.06s/it]
Evaluating coherence with PAI_R-FDR_longcalib


 57%|█████▋    | 4/7 [00:28<00:21,  7.16s/it]
Evaluating coherence with OIE_R-FDR_longcalib


 71%|███████▏  | 5/7 [00:35<00:14,  7.19s/it]
Evaluating coherence with BHM_R-FDR_longcalib


 86%|████████▌ | 6/7 [00:43<00:07,  7.25s/it]
Evaluating coherence with LMR_R-FDR_longcalib


100%|██████████| 7/7 [00:50<00:00,  7.21s/it]

We can now plot the results for each reconstruction. cohors indicate the level of coherency, while arrows show the phase angle (right = 0°, left = 180°, top = 90°, bottom = 270°, with intermediate directions representing angles in between).

period_ticks = np.logspace(1, 11, num=11, base=2)
for exp in exp_names:
    fig, ax = coh[exp].plot(yticks=period_ticks)
    ax.set_title(exp, fontweight='bold')
    
pyleo.closefig()
<Figure size 1000x800 with 2 Axes>
<Figure size 1000x800 with 2 Axes>
<Figure size 1000x800 with 2 Axes>
<Figure size 1000x800 with 2 Axes>
<Figure size 1000x800 with 2 Axes>
<Figure size 1000x800 with 2 Axes>

We see very little coherence aside from isolated “islands” that rise above the red noise benchmark, with phase angles (arrows) darting in all directions, suggesting no robust phase relationship. In some cases, there is a hint of an in-phase relationship (high coherency and phase angles close to 0°) at periods longer than 500 years. However, this signal lies near the cone of influence (the region of the scalogram affected by edge effects) and does not appear significant when tested against a red-noise (AR(1)) benchmark. The main exception is the BHM reconstruction from Barboza et al. (2019), which explicitly includes solar forcing as a predictor. It is therefore unsurprising that this reconstruction should exhibit coherence with solar variability, although the original study demonstrates that the signal is weak.

Phase Statistics

To further quantify this, we can use Pyleoclim’s phase_stats method, which characterizes the consistency of the phase relationship within a given scale band using two statistics from directional (circular) statistics:

  • σ (sigma): the circular standard deviation of the phase angles. A small σ indicates consistent phase angles across time; a large one suggests erratic phase behavior. Significance is assessed against an AR(1) null: if σ < σ_lo, the relationship is more phase-consistent than expected by chance.

  • κ (kappa): the concentration parameter of the Von Mises distribution — the circular analog of the normal distribution. Large κ means tightly clustered angles (strong phase lock). If κ > κ_hi, the phase relationship is significantly more concentrated than expected by chance.

We evaluate phase statistics over the centennial-to-millennial band (scales 100–2000 years), where solar forcing is thought to be most impactful.

rows = []
for exp in exp_names:
    ps = coh[exp].phase_stats(scales=[100, 2000], number=1000)
    rows.append({
        'Reconstruction'  : exp,
        'mean angle (°)'  : np.round(ps.mean_angle / np.pi * 180, 1),
        'σ'               : np.round(ps.sigma, 3),
        'σ_lo (5%)'       : np.round(ps.sigma_lo, 3),
        'σ < σ_lo?'       : ps.sigma < ps.sigma_lo,
        'κ'               : np.round(ps.kappa, 3),
        'κ_hi (95%)'      : np.round(ps.kappa_hi, 3),
        'κ > κ_hi?'       : ps.kappa > ps.kappa_hi,
    })

import pandas as pd
phase_df = pd.DataFrame(rows).set_index('Reconstruction')
phase_df
Loading...

The table above summarizes phase angle statistics for each GMST reconstruction vs. solar irradiance over the 100–2000 year band. Columns to focus on:

  • mean angle (°): the average phase lag; values near 0° indicate in-phase behavior (solar leads or lags by a small fraction of the period), values near ±180° indicate anti-phase.

  • σ<σlo\sigma < \sigma_\mathrm{lo}? and κ>κhi?\kappa > \kappa_\mathrm{hi}?: Boolean flags indicating whether the phase relationship is significantly more consistent than expected from a red-noise (AR(1)) null. Both flags must be True for the phase lock to be considered statistically significant.

Neither flag is satisfied for any reconstruction, confirming the absence of a robust, consistent phase lock between solar irradiance and GMST at centennial-to-millennial scales. Note, however, that spectral relationships are neither necessary not sufficient to establish causality. For a look at this question over a slightly longer timescale, the reader is referred to Landers et al. (2025).

References
  1. Neukom, R., Barboza, L. A., Erb, M. P., Shi, F., Emile-Geay, J., Evans, M. N., Franke, J., Kaufman, D. S., Lücke, L., Rehfeld, K., Schurer, A., Zhu, F., Brönnimann, S., Hakim, G. J., Henley, B. J., Ljungqvist, F. C., McKay, N., Valler, V., & von Gunten, L. (2019). Consistent multidecadal variability in global temperature reconstructions and simulations over the Common Era. Nature Geoscience, 12(8), 643–649. 10.1038/s41561-019-0400-0
  2. Barboza, L. A., Emile-Geay, J., Li, B., & He, W. (2019). Efficient Reconstructions of Common Era Climate via Integrated Nested Laplace Approximations. Journal of Agricultural, Biological and Environmental Statistics, 24(3), 535–554. 10.1007/s13253-019-00372-4
  3. Emile-Geay, J., McKay, N. P., Kaufman, D. S., von Gunten, L., Wang, J., Anchukaitis, K. J., Abram, N. J., Addison, J. A., Curran, M. A. J., Evans, M. N., Henley, B. J., Hao, Z., Martrat, B., McGregor, H. V., Neukom, R., Pederson, G. T., Stenni, B., Thirumalai, K., Werner, J. P., … Zinke, J. (2017). A global multiproxy database for temperature reconstructions of the Common Era. Scientific Data, 4(1). 10.1038/sdata.2017.88
  4. Landers, J. P., Emile-Geay, J., James, A. K., Munch, S. B., Khider, D., & Bard, E. (2025). A causal examination of the solar influence on Holocene climate. Wiley. 10.22541/essoar.176598901.17168843/v1
Frequency-based investigations of Paleoclimate Records
Frequency-based investigations of Paleoclimate Records
Exploring spectral properties of Holocene climate records
Querying the LiPDGraph