Best Pages 2k Proxies

By Andy May

Over the course of the past two posts, I have been searching through the original 647 Pages 2k proxies searching for the best proxies. As long-time readers will know, I am not a fan of statistically combining proxies into hemispheric or global temperature reconstructions. Proxies that adequately match historical, geological, and archeological records in the region around the proxy can be used to study local climate changes (Soon & Baliunas, 2003) & (Soon et al., 2003b) but not global or hemispheric changes. Region to region variation is too large, the few long-term proxies are too inaccurate, and the weather processes around the world are too poorly understood.

Due to loosely interconnected regional weather oscillations, such as the AMO (~64 years) and ENSO (~2.4 – 5.5 years), weather (and climate) varies both temporally and by region. Geological records (proxies) suggest that many of these oscillations have followed a roughly regular pattern for thousands, if not millions of years. It is possible that these oscillations stem from Earth’s orbital cycles combined with the sun’s variability and the orbits of the Moon and other planets in the solar system (Scafetta, 2010), (Stefani, 2026), (Stefani, 2021). How these forces interact to affect our climate is not fully understood at this time, but we can observe climate oscillations through well-located weather stations and proxies. There are two very high-quality proxy aggregates (a local reconstruction from multiple proxies) that I often refer to. They are an Indonesian ocean proxy by Yair Rosenthal (Rosenthal et al., 2013) and a Greenland proxy by Bo Vinther (Vinther et al., 2009). The two proxy aggregates are plotted in figure 1. Since most proxies are aggregates, I will use “proxy” and “proxy aggregate” interchangeably in this post, as is commonly done.

This post, combined with my last two posts (here and here) describe my search for more high-quality proxies that adequately match well-documented historical, glacial, and archeological records. A proxy’s credibility hangs on how well it matches the known history of the area and the glacial retreats and advances in the region. Regardless of its correlation with the modern instrumental record, if it does not match these local events, it must be considered invalid. Proxy temperature estimates and proxy chronologies often contain substantial uncertainties; thus, they need supporting outside data to be believed. Historical, archeological, and glacial records are useful benchmarks for judging proxy quality.

The defining features of figure 1 are the Holocene Climatic Optimum, from roughly 8000 BC to approximately 4000 BC, followed by Neoglacial cooling into the depths of the Little Ice Age. The coldest point in the Little Ice Age varies regionally, in Europe it is between ~1650 and 1715, but in the southern Pacific it is between ~1800 and 1820. Following the Little Ice Age, most of the world warms into the Modern Warm Period. The exceptions are most of the Southern Ocean, parts of the southern Pacific, and nearly all of Antarctica.

Some important historical events are not identified in figure 1, but the historical Bronze Age Warm Period (aka the Minoan or Mycenaean Warm Period) which brought the world the Egyptian New Kingdom, the Shang Chinese dynasty, and the Minoan Palace civilization, is clearly seen from about 1800 BC to 1000 BC with different peaks in each record. Likewise, the historical Bronze Age collapse (Cline, 2014), that occurs around 1200 to 800 BC can be seen in the records, followed by the Greek Dark Age, the Roman Warm Period,^[1] the European Dark Age (~400 to 800AD), the Medieval Warm Period (~800 to 1200AD), and the Little Ice Age. None of these warm and cold periods are globally synchronous, they generally vary in time by latitude and region and are most prominent in the middle to high northern latitudes. In general, they are more subdued in the tropics and Southern Hemisphere.

What was learned

In my last post, I listed what I learned from my initial inspection of the Pages 2k proxies, the list is summarized here:

1. Instrumental temperature/proxy regressions can perform poorly when the calibration interval is short, noisy, geographically mismatched, or characterized by weak proxy sensitivity. “Standard conversions” like Dansgaard et al. created for δ18O proxies or the Mg/Ca conversion by Anand & Elderfield are preferred.
   
2. Long-term climate changes can be obscured by regional variability. Regions change at different times and different rates, thus averaging them can remove major climate changes, like the Medieval Warm Period or the Little Ice Age.
   
3. The Medieval Warm Period and the Little Ice Age climate anomalies are much stronger in the Northern Hemisphere.
   
4. Combining proxies requires careful attention to regional coherence, historical events, and proxy selection because significant regional climate events may be attenuated or shifted in time. Good proxies are useful locally if they are supported by historical, geological, and archeological data.
   

Keeping these points in mind, I made one more pass through the remaining proxies. Figure 2 shows the location of the 255 non-tree-ring proxies as gray dots. The 128 “final” proxies are overlain in blue. The 57 “best” proxies are circled in blue and the subject of this post.

Proxy Selection Criteria

The most significant and well-documented climatic events in the historical record in the past 2,000 years are the Medieval Warm Period and the Little Ice Age (Behringer, 2010) & (Parker, 2012). These events (see this timeline) are very well documented in the Northern Hemisphere. Geological records (mainly fossil glacial moraines) show these events can also be seen in the Southern Hemisphere and tropics but not as clearly as in the Northern Hemisphere (Bray, 1968) & (Solomina et al., 2015). Neither event was synchronous around the world, so local or regional timing around each proxy record is important. In this post, the local timing and significance of each historical or geological event, if available, must be replicated by a proxy for it to be selected as one of the best proxies.

Next, I excluded all proxies that had fewer than four values between 1870 and 1900^[2] since this was my anomaly mean reference period and the earliest period with reasonable instrumental data. I also eliminated proxies that did not have valid measurement units. As summarized in Table 1, these restrictions removed an additional 127 proxies, leaving 128. Of the 128, I removed an additional 71 additional proxies for the reasons described in this post.

Stage	Remaining proxies
Page 2k Original	647
Remove tree rings	255
Remove poor coverage	128
Best proxies	57

Using “standard” correlations to convert proxies to ΔT

As discussed in the first post, many ice core researchers use standard correlation slopes (δ18O/°C) to convert their ice-core proxy δ18O values to temperature. These correlations are produced from temperatures and δ18O (or another proxy measurement) measured from the same sample, rather than comparing an air temperature or ocean temperature record to a possibly unrelated proxy. These sample-based slopes are valid today and vary regionally as shown in figure 3. Guan et al. note that while these slopes work today, they also vary temporally. Thus, they are not constant over time, and lack stationarity (Guan et al., 2016). This simply adds to the complexity and inaccuracy of calculating temperature from ice-core and other δ18O measurements. In this work I did not apply a temporal correction to the slopes and assumed stationarity, even though it isn’t so. This assumption is also made by Pages 2k and most researchers.

Correlations of proxies to unrelated instrumental temperatures often produce anomalous results. Examples of “standard” conversion factors by Dansgaard and colleagues (Dansgaard et al. 1964 & 1973) are shown in figure 3. It is important to consider that outside the areas shown in figure 3, current δ18O slopes can vary widely (Guan et al., 2016) and can often be negative in the tropics and other regions. This is on top of possible temporal problems with these slopes.

The present-day conversion factors given in figure 3 all have acceptable statistics, except for the North Atlantic locations, and are preferred to those derived from correlations. Pages 2k does not do conversions from δ18O, and other proxy measurements, to temperature individually. They use a complicated multiple proxy process that involves correlations between proxies and temperature over a wide area with a radius of 2,000 km. However, their methodology is still fundamentally a linear proxy to temperature correlation (Pages 2K, 2017) that requires a minimum of 20 proxy measurements over their instrumental period (1850-2000) and p <0.05.

My interest is in finding the highest quality individual proxies, so I did not attempt to duplicate their multiple proxy technique. Instead, I correlated each proxy with the points it had within the NOAA GlobalTemp instrumental grid (from 1870 to 2020) at the grid cell closest to the proxy. The OLS (ordinary least squares) correlations to the Vinther Agassiz and Renland ice core proxies are shown in figure 4 as the red curves in the upper plots. The blue curves are the ΔT conversion using the Coastal Greenland conversion factor from figure 3.

In figure 4 the upper plots show that the figure 3 conversion (blue) compares reasonably well to the OLS linear correlation (red) with instrumental data for the Renland ice core, but very poorly to the 2008 Agassiz proxy aggregate. Both records are straight from the Pages 2k database and their Pages 2k proxy name is in the title. These two locations are the principal raw data for Vinther’s 1979 aggregate Greenland proxy (Vinther et al., 2009) shown in figure 1 in orange. The lower plots in figure 4 compare the two ΔT conversions for each proxy. The ΔT from the OLS correlation is on the y axis and the standard ΔT from figure 3 is on the x axis. While the Agassiz correlation has better statistics, the visual correlation in the upper plot is poor, thus the conversion factor derived from OLS is likely flawed.

Comparing figures 4 and 5 helps explain the problem. The correlation-based conversion factor has a very low slope (0.044 °C/‰) as shown in figure 5. The OLS statistics can be misleading when the slope is low. Notice the large y axis scale difference between the two lower plots in figure 4. A correlation coefficient alone cannot determine whether a calibration is physically useful. A calibration with a very small slope may reproduce instrumental temperatures statistically while producing almost no meaningful proxy variability after conversion.

In figure 5, the low slope from the Agassiz conversion calibration plot (upper left) means almost any δ18O can match a NOAA temperature using this OLS fit. As a result, the converted ΔT in the upper right plot is nearly a straight line. The Renland slope (lower left) is steeper (1.1 °C/‰) and this results in a more active computed ΔT in the lower right plot that compares better with the standard conversion as shown in the upper right plot in figure 4.

Long-term climate changes can be obscured by regional variability

Depending upon the proxies chosen and how they are processed, major climate features such as the Medieval Warm Period, can appear or disappear in a composite, as shown in figure 6.

Figure 6 shows how sensitive major climatic features are to proxy selection, even when only the best proxies are used. This “averaging out” of major climatic changes is generally caused by events occurring at different times in different places due to regional climate oscillations, but poor dating of proxy samples can contribute to the averaging problem.

The Medieval Warm Period and the Little Ice Age climate anomalies

The Medieval Warm Period and the Little Ice Age are well-documented in historical accounts (Behringer, 2010) & (Parker, 2012). These climate anomalies are much stronger in the Northern Hemisphere. The proxies in figure 7 are just standardized temperature conversions, either using the Dansgaard temperature conversion slopes from figure 3 or the conversion used by the original proxy researcher (named “as is”).

In figure 7, using just the proxies converted using standard conversions or proxies converted to temperature by the original proxy researcher, the Medieval Warm Period and Little Ice Age are very prominent in the Arctic. The periods are noticeable, with considerable scatter, in the middle northern latitudes, and from the equator to 30N. The periods are very hard to see in the tropical southern latitudes, but visible in the middle southern latitudes.

For the plots in figure 7, gaps of up to 20 years are interpolated in each proxy to increase the chance that any given year has multiple proxies to average. Twenty years was chosen because the normal definition of a climatic event is one that lasts longer than 30 years. The means in figures 6, 7 & 8 are not area weighted since they are already segregated by latitude. They are also not weighted by uncertainty, as uncertainty is difficult to compute when regional variability is expected to be high.

All suitable proxies are plotted in figure 8, they were derived using OLS fits. With these proxies the Medieval Warm Period and Little Ice Age are visible in all latitude bands except for the Antarctic band from 90S to 60S, although both are subtle in the tropics.

Both figure 7 and 8 only show a very prominent Medieval Warm Period and Little Ice Age north of 30N. The latitudes north of 30N also display a lot of proxy scatter, suggesting a lot of regional variability.

Statistically combining proxies can attenuate climatic events

Figures 7 and 8 emphasize how risky combining proxies is. Proxy selection, even among the very best proxies can make a big difference in a statistical composite. For this reason, I still recommend using proxies to evaluate climate change one at a time. The highest quality proxies, in my opinion, are the proxies converted via Dansgaard’s conversion factors in figure 3 or converted to temperature by the original proxy researcher. Unfortunately, many proxies fall outside these groups.

Most of the proxies that did not make it into the “best proxy” list did so because their conversion from raw proxy units to temperature were judged to be unreliable. This was usually because the conversion slopes for that type (especially δ18O) were very inconsistent from proxy to proxy, the slopes were very small, or they did not match historical records in the area around the proxy.

Use the link below this post to download plots of the best proxies. Examples from the downloadable pdf are shown as figure 9 annotated with important historical events.

The upper left graph in figure 9 is from data collected and analyzed by Oppo, et al. (Oppo et al., 2009) from the Makassar Strait in Indonesia. The most severe drought in Indonesia’s history occurred between 1790 and 1810, the persistence of this drought well into the 19^th century was partially caused by volcanic eruptions from 1809 to 1835 in the area. The drought was accompanied by a series of very strong El Niño events and several monsoon failures (Rodysill et al., 2013).

The upper right graph in figure 9 is from a cave in Switzerland and data collected and analyzed by Larocque-Tobler, et al. (Larocque-Tobler et al., 2010). It shows the iconic European Medieval Warm Period in the beginning followed by the Little Ice Age. The profile is exactly as described by Behringer (Behringer, 2010) and Geoffery Parker (Parker, 2012), who provide considerable historical evidence for the temperature profile shown.

The lower left plot is from Antarctic ice core data collected and analyzed by Steig, et al. (Steig et al., 2000). This is one of the more active Antarctic proxies but still shows little variation over the past 2,000 years.

The plot on the lower right of figure 9 shows data from a 350-year-old stalagmite record in a South African cave collected and analyzed by Sundqvist et al., 2013. The δ18O record from this stalagmite correlates well with local temperature records and suggests a sharp drop in temperature in the early 1700s contemporaneous with a cold and dry period in South Africa (Hannaford & Beck, 2021). Dry spells in South Africa are often contemporaneous with El Niño events (Usman & Reason, 2004), which were relatively common in the early 1700s.

Glacial evidence

Glacial advances and retreats are “highly reliable indications” of climate change (Solomina et al., 2008). This is well documented in the literature (Bray, 1968) & (Solomina et al., 2015). The second-best climatic indicators are historical records (Soon et al., 2003b), (Behringer, 2010), and (Parker, 2012). The name for the long period of declining temperatures and advancing glaciers from the end of the Holocene Climatic Optimum around 4000BC into the depths of the Little Ice Age is the “Neoglacial” for a reason. To orient the reader, figure 10 shows the number of glacial advances, as well as the solar grand minima (SGM) and maxima (SGMx) for the past 2,000 years. This is qualitative data, but very reliable. The duration of the SGMs and the SGMxs is shown on the right-hand scale and is inverted so the SGMs have a shorter duration upwards, indicating colder at the bottom of the plot. Likewise, the number of worldwide glacier advances is inverted on the left-hand scale to show colder conditions downward.

Figure 10 shows that the longest Solar Grand Minimum (SGM) falls in Little Ice Age at about 1470 AD. This is also the time of the greatest number of global glacial advances. The maximum number of glacial advances globally (13) occurred from 1450-1550, with eight in the Northern Hemisphere, two in the tropics, and three in the Southern Hemisphere. Twelve occurred from 1550 to 1700, also more than any other period in Solomina’s tabulation (Solomina et al., 2015). The previous record was 10 global glacier advances during the Bronze Age Warm Period collapse around 1200 BC that marks the end of the Mycenaean Palace civilization.

Discussion

Bottom line, it is hard to rely only on geological temperature proxies. The data can be fairly accurate, but converting the proxy measurements to temperature is problematic. As discussed, the “standard” δ18O conversions are accurate today in their region but vary a lot from region to region and temporally (Guan et al., 2016). The proxies from Pages 2k are not very accurate, they need support from other sources, like historical records or glacial movements, and should probably not be combined except over small regions, that is converted into proxy aggregates.

Regional climate oscillations are not synchronous and often move in opposite directions. Regional climate changes complicate global averages. This also means that major climate changes, like the Little Ice Age and the Medieval Warm Period, while very visible in regional temperature reconstructions can be “averaged out” globally or hemispherically. Important qualitative historical and glacial records should not be ignored; these provide support or destroy confidence in proxies.

Download plots of the best proxies here.

Download the R code here.

Download the bibliography here.

Appendix: The R code used for this analysis

This R code zip file can be unzipped in the same working directory that contains the contents of the R code zip file from my last post. The “root_dir” variable in these R programs must be set to the same working directory on your computer. The programs discussed in this appendix are:

D18O_standard.r (Makes most of the figures)

D18O_standard_slope.r (makes figure 3)

Map_proxies.r (makes figure 2)

The zip file also contains some necessary data and a revised version of Pages2k.r, described in the previous post. The changes to Pages2k.r are minor and will not be described, they are mainly to create final_proxies.RData, which is used by δ18O_standard.r.

d18O_standard.r

This program was designed to compare “standard” δ18O temperature conversions (Dansgaard et al., 1964, 1973) to conversions from OLS fits to NOAA modern instrumental data. It also looks at overall conversions for ocean δ18O data. Ocean δ18O values from sediments and corals are often used to estimate paleotemperatures, but the slopes are all over the place, both positive and negative, and not very useful. I failed to find a good overall ocean conversion factor and eventually abandoned the idea, but left the code in the program in case others want to pursue this problem. Until someone figures out a standard way to convert δ18O to temperature in the ocean, use the original proxy researcher’s chosen method.

The program first loads final_proxies.RData. It contains the 128 final proxies discussed in the previous two posts. Then after the long code block examining and plotting the ocean proxies, the program creates a pdf file “Proxy_DeltaT_Comparison.pdf” and populates it with plots that compare the standard δ18O conversions to temperature to the OLS best fit conversions for each of the best 57 proxies. Right after creating two lists, “delta_T_S_results” (holds the standard delta_T) and “best_proxies” (the list of 57 best proxies), the pdf file is created with the following code:

# Create pdf
pdf(
  paste0(plots_dir,"Proxy_DeltaT_Comparison.pdf"),
  width = 9,
  height = 5
)

Then the program enters a loop through the proxies. All proxies are then turned into anomalies from the mean of 1870 to 1900:

###########################################################
# CONVERT TO ANOMALIES
###########################################################

  delta_T_S <- make_anomaly(
    year,
    T_std,
    1870,
    1900
  )

Then the plots of original and standard delta_T, see figure 4, are made. They are written to the pdf file.

plot(
  year,
  proxy$delta_T,
  type = "l",
  lwd = 2,
  col = "red",
  ylim = common_range,
  xlab = "Year",
  ylab = expression(Delta*T~"(°C)"),
  main = paste0(
  proxy_name,
    "\nMethod: ",
    method, ", Type: ", Type
  )
)

lines(
  year,
  delta_T_S,
  col = "blue",
  lwd = 2
)

abline(h = 0,
  lty = 2,
  col = "gray50")

legend(
  "topleft",
  legend = c(
   "Original delta_T",
   "Standard delta_T_S"
  ),
  col = c("red", "blue"),
    lwd = 2,
    bty = "n"
)

# Fit OLS model
fit <- lm(proxy$delta_T ~ delta_T_S)

# Extract statistics
slope <- coef(fit)[2]
r2    <- summary(fit)$r.squared
  best_proxies[[length(best_proxies)+1]] <-
    data.frame (
      site = proxy$site,
      latitude = proxy$latitude,
      longitude = proxy$longitude,
      elevation = proxy$elevation,
      Region = proxy$Region,
      meas_type = proxy$meas_type,
      meas_units = proxy$meas_units,
      meas_archiveType = proxy$meas_archiveType,
      slope = slope,
      R2 = r2
    )


# Create plot
plot(
  delta_T_S,
  proxy$delta_T,
  type = "p",
  pch = 4,
  cex = 0.7,
  col = "red",
  xlab = "Standard delta_T_S",
  ylab = "Original delta_T",
  main = paste0(
    proxy_name,
    "\nMethod: ",
    method, ", Type: ", Type
  )
)

# Add OLS line
abline(fit, col = "blue", lwd = 2)

# Add annotation
legend(
  "topleft",
  legend = c(
    sprintf("Slope = %.3f", slope),
    sprintf("R² = %.3f", r2)
  ),
  bty = "n"
)

Next a list of the best proxies is read, note this list can be changed:

# ----------------------------------------------------------
# Read preferred proxy list
# ----------------------------------------------------------

best_proxy_list <- read.csv(
  "csv/Best_proxies_final_list.csv",
  stringsAsFactors = FALSE
)

use_sites <- unique(best_proxy_list$site)

And gaps in each proxy that are less than 20 years are interpolated. Twenty years is chosen because it is the longest period that will not affect the detection of climatic effects (~30 years):

for(site in intersect(use_sites, names(delta_T_S_results))) {
  
  x <- delta_T_S_results[[site]]
  
  df <- data.frame(
    year = round(x$year),
    delta_T_S = x$delta_T_S
  )
  
  # One value per year
  df <- df %>%
    group_by(year) %>%
    summarise(
      delta_T_S = mean(delta_T_S, na.rm = TRUE),
      .groups = "drop"
    )
  
  yrs_full <- seq(
    min(df$year, na.rm = TRUE),
    max(df$year, na.rm = TRUE),
    by = 1
  )
  
  interp_df <- merge(
    data.frame(year = yrs_full),
    df,
    by = "year",
    all.x = TRUE
  )
  
  interp_df$delta_T_S_interp <- na.approx(
    interp_df$delta_T_S,
    x = interp_df$year,
    na.rm = FALSE,
    maxgap = 20
  )
  
  delta_T_S_results[[site]]$delta_T_S_interp <-
    interp_df$delta_T_S_interp
  
  delta_T_S_results[[site]]$interp_year <-
    interp_df$year
}

Then the latitude band plots are made. This time, due to interpolation, the chances that multiple proxies have values in a given year are improved. In this code the plots seen in figure 7 are made, among others.

p <- ggplot() +

# --------------------------------------------------------
  # Plot
  # --------------------------------------------------------
  
  p <- ggplot() +
    
    geom_line(
      data = proxy_df,
      aes(
        x = year,
        y = delta_T_S,
        group = site
      ),
      color = "grey75",
      linewidth = 0.3,
      alpha = 0.8
    ) +
    
    geom_ribbon(
      data = summary_df,
      aes(
        x = year,
        ymin = lower,
        ymax = upper
      ),
      fill = "lightblue",
      alpha = 0.4
    ) +
    
    geom_line(
      data = summary_df,
      aes(
        x = year,
        y = mean_deltaT
      ),
      color = "blue",
      linewidth = 1.2
    ) +
    
    coord_cartesian(ylim = c(-3.0, 3.0), 
                    xlim = c(0, 2000)) +
    
    labs(
      title = paste(
        "Standardized Proxy Composite:",
        gsub("_", " ", band_name)
      ),
      subtitle = paste0(
        length(unique(proxy_df$site)),
        " proxies"
      ),
      x = "Year",
      y = expression(Delta*T[S]~"(degC)")
    ) +
    
    theme_bw(base_size = 18)
  
  outfile <- paste0(
    plots_dir,
    "Latitude_Band_delta_T_S_",
    band_name,
    ".png"
  )
  
  print(p)
  
  ggsave(
    outfile,
    p,
    width = 12,
    height = 7,
    dpi = 300
  )
  
  cat("Created:", outfile, "\n")

Δ18O_standard_slope.r

This program reads the digitized Dansgaard δ18O data computes OLS lines and statistics for each region and plots them. One of the plots is shown in figure 3. It also produces a regression summary table.

Map_proxies.r

This program produces the map shown in figure 2. The map uses the Robinson projection. A latitude/longitude grid (or “graticule”) is created and plotted, and labels for both longitude and latitude are produced. Each set of points (all non-tree-ring proxies, final proxies, and best proxies) are plotted separately.

1. Bianchi & McCave, 1999 give the peak time at 100 AD but dates vary in historical accounts and in proxies. If there is one peak, it is probably between ~400 BC to ~100 AD or there may be two peak Roman Warm Period temperatures most places, hard to say. The warm period exists, but it is variable by location. ↑
   

   
   

2. Prior to 1870 the uncertainty in SST (sea surface temperature) is simply too high to use and oceans cover 70% of Earth’s surface (Kennedy et al., 2011) & (Kennedy et al., 2011b). ↑