The Pages 2k R code and Temperature Records by Latitude Band

By Andy May

The description of the Pages2k.R program is in the Appendix at the end of this post.

In my last post, I discussed the problems with the Pages 2k temperature reconstruction since 0 AD using the Vinther Greenland ice-core proxies (Vinther et al., 2009) as an example. However, many of you may want to inspect more of the Pages 2k proxies, so I am making my R code available for download and provide some documentation of critical parts of it in the appendix to this post. The R code and some of the data necessary can be downloaded here. To get the most out of the appendix and this post, I recommend downloading the code so you can view it as you read.

Making Latitude-band Composite Plots

Most of the Pages2k.R code just examines the non-tree-ring proxies in the database and creates basic quality control plots, however climate varies significantly by latitude (Scotese et al., 2021). Due to the seasons and Earth’s orbit around the sun, each latitude band can change its climate at a different rate and not always in the same direction. So, after the individual proxies are processed and plotted, the last section of the R code in Pages2k.R processes the proxies in latitude bands. The program divides the globe as follows:

# Latitude bands
lat_bands <- list(
"90N_to_0" = c(0, 90),
"90S_to_60S" = c(-90, -60),
"60S_to_30S" = c(-60, -30),
"30S_to_0" = c(-30, 0),
"0_to_30N" = c( 0, 30),
"30N_to_60N" = c( 30, 60),
"60N_to_90N" = c( 60, 90)
)

Pages2k has a for loop to collect all the final proxy series into each band and plot them. These are all the suitable non-tree-ring proxies in the database. Figure 1 shows four of the plots:

Figure 1. Plots of all the proxies in final_proxies for 60N to 90N, 90S to 60S, 60S to 30S, and 30S to the equator. Tree ring proxies are excluded as well as proxies that were excluded by Pages 2k, or had insufficient data in the calibration period from 1870-2020.

The proxies plotted in figure 1 have all been converted to temperature anomalies from the 1870-1900 mean. The proxies in each plot are searched for years that contain 3 proxies or more. If encountered the mean and standard deviation (sd) are computed as well as the mean + sd and mean – sd. The mean is plotted in blue and the one standard deviation region is shown with light blue shading. Following is the code that does all this:

summary_df <- proxy_df %>%
group_by(year) %>%
summarise(
mean_deltaT = mean(delta_T, na.rm = TRUE),
sd_deltaT = sd(delta_T, na.rm = TRUE),
n = sum(!is.na(delta_T)),
.groups = "drop"
) %>%
mutate(
mean_deltaT = ifelse(n >= 3, mean_deltaT, NA),
sd_deltaT = ifelse(n >= 3, sd_deltaT, NA),
lower = mean_deltaT - sd_deltaT,
upper = mean_deltaT + sd_deltaT
)

The use of the R programming pipe or “%>%” and “summarise” are discussed in this post. This R operator is very useful when processing climate data. For more details about the R code see the appendix.

Latitude band plot discussion

Figure 2 shows the remaining latitude bands. The plots are illuminating.

Figure 2. The remaining “final_proxies” latitude band plots. These are for the equator to 30N and 30N to 60N. As in figure 1, tree rings are excluded.

The number of proxies in each latitude band is variable, 60S to 30S only has three proxies and 30S to the equator has 34. The number of proxies and proxy aggregates in each band must be kept in mind when interpreting these plots. In the Arctic (60N-90N), the warmest period (judged only from those years with 3 or more proxies shown in blue) was from 750 AD to 850AD when the Vikings began their conquest of the British Islands and Ireland. Iceland was settled by the Vikings in 874 and Greenland in about 982. According to these limited proxies, this entire period was warmer than today in the Arctic. This is very limited data but it does make sense historically.

In contrast, Antarctic temperatures (90S-60S) were flat according to nearly all the proxies. The Southern Hemisphere (60S-30S) shows the Medieval Warm Period as warm as today where the three proxies agree. It shows a noticeable Little Ice Age from around 1400 to 1600 and second cool period around 1800. Yet the middle southern latitudes have less dramatic changes than the northern middle latitudes in figure 2.

The southern tropics from the equator to 30S also show a modest Medieval Warm Period and a subdued Little Ice Age but the Little Ice Age has moved later in time and corresponds better with the Northern Hemisphere.

The band from the equator to 30N in figure 2 is a bit of a mess, but the longer-term proxies show the Medieval Warm Period to be warmer than today, and this is certainly true for the mid Northern Hemisphere from 30N to 60N in figure 2 on the right.

Longer-term proxies only

The next set of plots only includes longer-term proxies with data from 200 AD or before and from 1900 AD or later. These are shown in figures 3 and 5. This proxy subset shows a flat trend in both the Antarctic and the Arctic, but for different reasons. In the Antarctic all the proxies agree that the trend is flat, except for the borehole temperature survey shown with the long, smooth, loopy record. In the Arctic the individual proxies are very wild and the flat trend in the mean appears to be due to opposing proxies “averaging out.”

In the background, the individual Arctic proxies show the familiar Medieval Warm Period and cooler Little Ice Age. This limited data shows the Arctic is much more active climatically than the Antarctic. The middle Southern Hemisphere band is only represented by one proxy which is anomalously warm in the Little Ice Age and anomalously cool in the Medieval Warm Period. However, the equator to 30S band shows both the Medieval Warm Period and the Little Ice Age in their respective places, they are just a little more subtle than in the Northern Hemisphere.

Figure 3. Only records with data from 200 AD or before and to 1900 AD or later.

The lone proxy in 60S-30S is “SAm-LagunaChepical.deJong.2013.” These are late summer temperatures (Nov to Feb) from a high elevation lake in a sub-tropical (32°S) Andes glacial lake in Chile. The record was constructed from visible spectrum reflectance spectroscopy, and the record reasonably correlates to NOAA GlobalTemp (R2=0.24 & p=1.7 x 10-9). However, the lake is over 3,000 meters in elevation, it is a summer proxy, and a glacial lake, thus it is not well located as a proxy for the entire Southern Hemisphere middle latitudes. Further, according to the Proxy 2k report summer proxies often do not reflect annual temperature trends in the middle and high latitudes.

The lone long-term proxy in 30S to the equator is “Ocn-MakassarStrait.Oppo.2009” (Oppo et al., 2009). The proxy has a distinct Medieval Warm Period and Little Ice Age. The Medieval Warm Period is slightly warmer than today in this proxy. This is an interesting choice of proxy; it is the earliest of three possible proxy aggregates of the same ten Indonesian sediment cores in the Makassar Strait area. The others are Rosenthal et al., 2013 and Linsley et al., 2010. Both Linsley et al. and Oppo et al. present SST (sea surface temperature) proxy aggregates that are complicated by varying flows through the Makassar Strait from the North Pacific, South China Sea, and Indian Oceans. All three proxies are compared in figure 4.

Figure 4. Three Makassar Strait proxy aggregates compared. The black line is the Rosenthal 500-meter depth temperature anomaly, shown with its 1 standard deviation bands as gray shading. Both Oppo and Linsley provide two SST aggregates, and both are shown. The red Linsley aggregate is nearly hidden by the green, as they are nearly identical. All these studies are by the same team and use the same set of sediment cores from the Makassar Strait area but are based on different foraminifera. All aggregates are anomalies from the 1700-1780 mean. Data sources: (Oppo et al., 2009), (Linsley et al., 2010), and (Rosenthal et al., 2013).

The aggregate proxies in figure 4 are all made from a standard formula relating Mg/Ca to temperature (Anand et al., 2003) and are not derived from a correlation to GlobalTemp or HadCRUT4. This is analogous to the Dansgaard (1964 & 1973) “standard” Greenland temperature conversion (Dansgaard, 1964) & (Dansgaard et al., 1973) mentioned in the previous post. Pages2k did not use the raw Mg/Ca data, they used the temperatures provided by Oppo et al. as shown in purple in figure 4 and shown in figure 3 (lower right). Notice that figure 4 is shown with the original 10-year samples from Oppo (2009) and the sample interval used by Pages 2k is variable, but about 5 years, this changes some of the details but clearly they are based on the same data. The Linsley reconstruction (green in figure 4) was also in the Pages 2k database, but was not selected since it stops in 1780 before the instrumental calibration period. Due to the Linsley data, the reference mean used in figure 7 has been moved to 1700 to 1780, which shifts the Oppo curve (in purple) upward about three quarters of a degree from figure 3.

This region is within the Indo-Pacific Warm Pool (IPWP), which has a current annual average SST of about 29°C, near the maximum possible open ocean surface temperature of about 30°C (Sud et al., 1999). SSTs in this region are an important factor in controling planetary-scale atmospheric circulation, and together with polar temperatures, exert an outsized influence on global climate (Oppo et al., 2009). Of the proxy aggregates provided by Rosenthal et al., is an aggregate of the ocean mixed layer (up to 500 meters depth) temperature in the Makassar Strait. He used foraminifera fossil shells that calcify deeper in the ocean than those used by Oppo and Linsley. At this depth the predominant flow is from the North Pacific Ocean to the Indian Ocean and it is representative of North Pacific SSTs on climatic time scales. The various aggregates from Oppo (2009) and Linsley (2010) are SST aggregates from surface or shallow dwelling formainifera in the Makassar Strait area, but as you can see in figure 4 they show a similar pattern to the deeper reconstruction by Rosenthal (2013).

Overall, Rosenthal’s aggregate proxy is probably better to use from a climatic perspective because it represents a large area in the North Pacific and sits in or near the IPWP which is a key component of long-term climate changes (Oppo et al., 2009). Figure 5 shows the remaining two long-record proxy collections.

Figure 5. Long-term Pages 2k proxy aggregates for the equator to 30°N and 30°N to 60°N. The blue bars in the right-hand graph are years with three or more proxy values.

The lone proxy in 0-30N is “Ocn-EasternTropicalNorthAtlantic.Kuhnert.2011.” This is an aggregate proxy of eastern tropical Atlantic SSTs (Kuhnert & Mulitza, 2011). The cores used in the study are from offshore Mauritania and strongly reflect (95% confidence) the AMO (Atlantic Multidecadal Oscillation) ~64-year oscillation. Overprinted on the AMO pattern, the SSTs in this very limited dataset show a drop from the Medieval Warm Period around 1000 AD into the Little Ice Age from ~1300 to 1850 before warming into the Modern Warm Period, which, at this location, is cooler than the Medieval Warm Period.

The five long records between 30N and 60N show a much larger drop in temperature from the Medieval Warm Period into the depths of the Little Ice Age. The more extreme records in the northern middle latitudes show more than a two degree drop. This is a very limited set of proxies, but it does match historical records.

Limitations of this Analysis

This analysis is exploratory; I wanted to find out what data was used by Pages 2k and look at it proxy-by-proxy. It was never intended to be a formal analysis of the proxies by latitude band or region-by-region. I wanted to see the nuts and bolts. As a result, I skipped over a lot, in particular:

  • Tree-ring proxies were excluded by choice.
  • Proxy calibration uses nearest-grid NOAA temperatures.
  • OLS regressions were used for calibration; there are other techniques that could be used.
  • Proxy sample age uncertainties are ignored, yet they probably exist.
  • Latitude-band composites are simple averages, with no consideration that sample increments between proxies may not match. No interpolation was done, which would help the analysis.
  • Some latitude bands contain very few long records. Long records are preferred since they are consistent from beginning to end and not spliced.

Conclusions

One overall conclusion from this limited data is the Antarctic proxy aggregate shows substantially less multi-centennial variability than the Arctic proxy aggregate over the past 2,000 years (see figures 1 & 3). This implies little climate change in the Antarctic over the period, whereas the Arctic (upper left, figures 1 & 3) and IPWP (figures 1 & 3, lower right) have undergone significant changes over the same period. It appears global climate is influenced more by the Northern Hemisphere and the tropics than the Antarctic. This is also suggested by the correlation between the detrended AMO and global mean surface temperature.

It is also possible that the middle Southern Hemispheric latitudes have experienced changes opposite to those seen in the Northern Hemisphere (compare figures 1 and 3). This pattern of the Southern Hemisphere (or portions of it) acting opposite to the Northern is also seen over the entire Holocene as shown in figure 6. This correlates with Earth’s orbital precession, which causes a different pole to face the sun during the summer roughly every 10,000 to 11,000 years at perihelion when Earth is closest to the sun.

Figure 6. Rosenthal’s Makassar Strait aggregate proxy, compare to the close up in figure 7, both in black, compared to the Antarctic Dome C ice core proxy and the Vinther 2009 proxy from the previous post.

The data in this post and the previous post suggest several things:

  1. Using correlations to calibrate a proxy to the instrumental era does not work well. Instead, more established standard conversions like Dansgaard et al. created for d18O for several regions or the Mg/Ca conversion by Anand & Elderfield are preferred. The instrumental period overlap with proxies is too short and the period can be anomalous.
  2. Long-term climate changes can be obscured by regional variability, as seen in the 60N to 90N plots in figures 1 and 3. Fluctuations, either caused by inaccurate dating or natural regional variability, can remove major long-term climatic changes. Even the Medieval Warm Period or the Little Ice Age can be removed by averaging.
  3. The Medieval Warm Period and the Little Ice Age do exist nearly everywhere, except in Antarctica. These long-term climate anomalies are strongest in the Northern Hemisphere and tropics but they are not globally synchronous.
  4. Combining proxies is very risky. Good proxies are relevant locally if they are supported by historical and archeological data, but they are not accurate or precise enough in either the temperature or time dimensions to be combined. Examining them one at a time with as much local historical data as possible is more illuminating, combining disparate proxies hides important details.

The R program described in the appendix below was mostly used to examine the Pages 2k database and make plots of the non-tree-ring proxies. The next post in this series will present several more focused R programs designed to reanalyze the Pages 2k proxies.

The R code can be downloaded here.

The bibliography can be downloaded here.

Appendix: Running the R code

If you are new to R, then download and install it from CRAN as well as RStudio. Then create an empty working directory on your hard drive and unzip my R code zip file (“Pages2k_R_code_Andy_May_v2.zip”) into it. After starting RStudio run “Install critical libraries.r,” it installs the R libraries that you need to run the Pages 2k programs.

The first R Program and downloading the data

After installing the R libraries, load “Pages2K.r” into RStudio. RStudio runs the statement at the cursor when you hold down the crtl key and press enter. You can also highlight multiple statements and press ctrl-enter to run all of them, but be careful, the selection will stay active after they run and pressing almost any key will delete the selected code. So immediately click in the blue highlighted code to clear the selection for safety.

After the program loads the necessary libraries, it sets “root_dir”, which is the working directory. Change the directory in the code to whatever directory contains the Pages2k R code. This is where R will do its work. A directory will be created to receive the plot output (plots/) and another to create text output (csv_files/).

Next Pages2k.r downloads NOAA’s GlobalTemp grid and the Pages2k RData proxy database. The GlobalTemp data is in the variable “nc.” The Pages2k database is stored in a variable called “D,” then Pages 2k checks its length with:

length(names(D))

The number printed should be 647, for 647 proxy files. Next Pages 2k makes an inventory of the proxy types and stores it in a data frame called “Inventory.” Below is the inventory you should see in the RStudio console window:

print(Inventory)
archiveType Count
1 Wood 392
2 Coral 82
3 MarineSediment 57
4 GlacierIce 47
5 LakeSediment 42
6 Documents 15
7 Sclerosponge 4
8 Speleothem 4
9 Borehole 3
10 Other 1

ArchiveType is a top-level variable in each site record in “D,” it categorizes the proxy. One of the sites is in central Chile, its Pages 2k name is: “SAm-LagunaAculeo.vonGunten.2009.” The original study is by Lucien von Gunten and colleagues (von_Gunten et al., 2009). It is an austral summer (DJF) temperature reconstruction based on reflection spectrometry data from a Laguna Aculeo sediment core. It has a robust correlation with modern instrumental records as shown in figure 8 below. The variable names for this site are:

names(D[["SAm-LagunaAculeo.vonGunten.2009"]])
[1] "paleoData" "agesPerKyr" "archiveType" "createdBy"
[5] "dataContributor" "datasetId" "dataSetName" "datasetVersion"
[9] "googleDataURL" "googleMetadataWorksheet" "googleSpreadSheetKey" "lipdverseLink"
[13] "maxYear" "minYear" "originalDataURL" "tagMD5"
[17] "changelog" "pub" "geo" "@context"
[21] "lipdVersion"

Individual commands can be typed into the console window of RStudio, the contents of archiveType can be inspected in this way:

D[["SAm-LagunaAculeo.vonGunten.2009"]]$archiveType
[1] "LakeSediment"

As shown in the inventory, there are a total of 647 proxies and 392 (61%) are “Wood” or tree ring proxies. As explained the last post, tree ring proxy series are detrended, so they are unreliable long-term climate indicators (Briffa et al., 1998b). This is the “segment length curse” of tree ring proxies (Cook et al., 1995). The Pages 2k team uses “Regional Curve Standardization” to try and regain some long-term signal from the tree rings, but since tree populations in any given area change over time this is an unreliable technique (Esper et al., 2003). For this reason, I chose to eliminate tree ring proxies from the database. You are quite free to change the program to retain them for your study if you wish.

After making the inventory, Pages2k.r downloads the NOAA NCDC version of the Vinther data that is used by Pages2k and places them in the subdirectory “NCDC.” This version is higher resolution (up to 1-year) than Vinther’ original proxies from iceandclimate.nbi.ku.dk. The original proxies, from the Denmark site, are pre-downloaded and in the “Proxies” subdirectory of root_dir.

Inspecting the Pages 2k contents

This is the first program I wrote in this study, so a lot of it is just inventorying the data and counting the number of proxies that fit certain criteria I thought were important when I started the study. This inspection process uses some repetitious code that could be shortened, but I left it in to help the reader do their own searches and inventories.

After downloading the NCDC proxies, the next statement is “names(D)” which lists all 647 proxies to the console.

As described above, to investigate the contents of any of them, set the “site” variable to the proxy name and place a “$” after “D[[site]]” and the variable name. “D[[site]]$pub” returns (among other publication details) for site =”Arc-Agassiz.Vinther.2008″:

D[[site]]$pub
[[1]]
[[1]]$author
[1] "Vinther, B. M.;Clausen, H. B.;Fisher,
D. A.;Koerner, R. M.;Johnsen, S. J.;Andersen, K. K.;
Dahl-Jensen, D.;Rasmussen, S. O.;Steffensen, J. P.;
Svensson, A. M."
[[1]]$doi
[1] "10.1029/2007JD009143"
[[1]]$journal
[1] "Journal of Geophysical Research"
[[1]]$pubYear
[1] "2008"
[[1]]$title
[1] "Synchronizing ice cores from the Renland and Agassiz
ice caps to the Greenland ice core chronology"

Which is very useful information, D[[site]]$geo$longitude and D[[site]]$geo$latitude return:

D[[site]]$geo$longitude
[1] -73.1
> D[[site]]$geo$latitude
[1] 80.7

Next notice the following statements:

pd <- D[[site]]$paleoData[[1]]
mt <- pd$measurementTable[[1]]
geo <- D[[site]]$geo[[1]]

These store part of the paths to variables in the critical paleodata, measurementTable, and geo sections of “D” to save typing. This is how mt$d18O$”useInGlobalTemperatureAnalysis”, which is shorthand for:

D[[site]]$paleoData[[1]]$measurementTable[[1]]$"d18O"$"useInGlobalTemperatureAnalysis"

is stated using these shortcuts. Thus:

mt$d18O$"useInGlobalTemperatureAnalysis"
[1] "TRUE"

The value “TRUE” means the proxy was used in the Pages 2k reconstruction. If the correlation statistics of the proxy to the modern instrumental era are in the database somewhere, I did not find them.

Functions

Two functions used in this program are listed next. The first is a function to make sure that the site names are safe to use as variable names in R (it eliminates certain dangerous characters from the name like “/” and “\\”). The second is a function to plot proxies called “plot_proxy.” Just press ctrl-enter to execute the functions, they will be syntax checked and stored for future use.

Initial processing

Following plot_proxy, a list called “rows” is created to store critical variables for each proxy. The data we loaded from Pages 2k, is from the list variable “D.” After declaring rows, Pages2k.r uses it in a “for” loop (“for (site in names(D))”) that loops through all the proxies in D and picks up key data items, like archiveType, latitude, longitude, and elevation. It also excludes all tree ring (archiveType = “wood”) proxies. It uses this data to build a list (“rows”) with the selected variables and then binds rows into a single data frame called proxy_df. A data frame is a special type of list where all components (aka columns) are vectors of the same length. At the end of the code block the data frame made from rows (proxy_df) is written out as a csv file to csv_dir in a file named “Proxies_without_tree_rings.csv”.

Examining “Proxies_without_tree_rings.csv” we see that it contains information on 255 proxies, including their latitude, longitude, elevation, type of proxy, region, and the title of the original paper describing the proxy.

The next code block also excludes tree rings but then searches the proxy “measurementTables” for the different types of measurements and counts them. The dominant type is d18O (113) and second is temperature (108).

This is followed by a code block that goes through the list of proxies again and rejects those that are not used in the Pages 2k reconstruction or have no proxy values in the measurement table (“mt”). If multiple variables are flagged as TRUE, it takes d18O or temperature. These checks do not change the number of proxies retained, the total is still 255, which is useful information. All this rejection logic can be changed; these choices are strictly opinion.

Uncertainty

The program does check for researcher supplied uncertainty. One of the proxies (“Afr-ColdAirCave.Sundqvist.2013”) had author supplied uncertainty, so this logic was included to capture others, but of the proxies I examined there was only one.

Year

Most of the proxies had time supplied as “bp” or years before 1950. If the proxy has time in these units Pages2k converts it to calendar_years.

Geographic information

The code block picks up latitude, longitude, elevation, region, and other details for each proxy.

Inventory

The next code block shows the number of remining proxies (“length(proxy_list)” which should still be 255, lists the first ten, and shows the number of rejected sites (0, except for the tree ring proxies). Then an inventory is built and written out to the csv_dir as “Proxy_list_Inventory”. The inventory is interesting, it shows that coral d18O, marine sediment, and lake sediment proxies make up more than half of the proxies that passed our filters.

Plots

The next code block makes plots of the surviving 255 proxies. Figure 7 is an example plot from plot_proxy. It is a rather distinctive plot of a South American temperature proxy.

Figure 7. This proxy is from central Chile, in the Southern Hemisphere. It is a summer (DJF) temperature proxy aggregate from reflection spectrometry. It has a robust correlation to summer temperatures in the instrumental record from 1901-2000 (von_Gunten et al., 2009).

In figure 7 it is interesting that the initial large temperature drop (at least in Chile) occurs just a few decades after the Wolf Solar Grand Minimum (SGM). Usoskin’s Little Ice Age Solar Grand Minima are located on the post (Usoskin et al., 2016). According to Usoskin, the Spörer SGM is the longest in over 8,000 years.

Computing ΔT from the proxies

After making the proxy plots and examining them, we begin the process of converting them into temperature anomalies. The first step is to build a data frame, called proxy_summary. It contains some critical variables like the site name, minimum year and maximum year, and a flag to the reference period of 1870-1900. While working with the proxies, I found some errors, these are fixed in the code block headed with the comment “Fix erroneous records, make changes.” After fixing the errors, the data frame final_proxies is created, it has a length of 128 proxies. All the proxies in final_proxies have valid units, the minimum year is 1880 or earlier and the maximum year is after 1900, and it has at least 4 values between 1870 and 1900.

In step 2 Pages 2k sets up to extract the NOAA GlobalTemp temperature record from the grid cell closest to the proxies with the functions nearest_grid and get_noaa_annual. This is used to calibrate the proxies in the calibration period 1870 to 2020. The NOAA GlobalTemp 5°x5° latitude/longitude grid cells are very large (1.2 million square kilometers at the equator) so over complex terrain or in ocean regions with strong ocean currents, surface temperature can vary a lot across one grid cell.

Next, in step 3 we enter a loop through all the final proxies. In the loop, we extract the NOAA GlobalTemp record for the site’s location and create a 5-point smoothed version of the proxy. Next the logic separates those proxies already in temperature units from those that are not. If the proxy is in temperature units, it was converted by the original researcher. In that case, we simply move it to an anomaly from the 1870 to 1900 mean and compare it to the NOAA GlobalTemp record and use it as is. The logic in doing this is that if the original researcher had a preferred way to convert his proxy to temperature, he is closest to the data and has already chosen the best methodology and temperature series to use. Our method is a crude simple regression against a 5×5 degree GlobalTemp grid. Mixing our crude method with the original researcher’s conversions is not ideal, although it is worth noting that the R2 values for researcher temperature vs NOAA GlobalTemp (R2 mean 0.19) are no better than for the original proxy and NOAA GlobalTemp (R2 mean of 0.15).

For the non-temperature proxies, the comparison to GlobalTemp is done with five-point smoothed data. After doing the comparison, we make a comparison plot, like in figure 8.

Figure 8. The instrumental era calibration plot for von Gunten’s spectral reflectance proxy in central Chile.

If the proxy is not already a temperature, but is in d18O or some other units, it is converted in Step 4. This is done using a linear transform after 5-point smoothing the data. Note the smoothing technique used here does not take the sample increment (or spacing) into account, so it varies temporally with sample spacing. The d18O to NOAA GlobalTemp calibration plot is shown for the Renland core in figure 9.

Figure 9. The instrumental era calibration plot for the Renland ice core, it has unusually good statistics.

Once the calibration is done, we move on to step 5 where the “predict” function uses the instrumental era fit statistics to compute temperature along the entire proxy record. Then the temperature is moved to an anomaly from the 1870-1900 mean. The calibration statistics are stored in final_proxies, and Pages 2k plots the resulting temperature record.

After all the sites are processed, another inventory is created and written out as “final_proxies_Inventory.csv.” This inventory reveals that coral d18O is the most common proxy (40) and glacier ice d18O is the second most common (31), followed by marine and lake sediment cores already converted to temperature (21).

Then Pages 2k counts the proxies by measurement type and reveals that 77 of the proxies are d18O measurements and 31 have been converted from their original units to temperature units already. This is followed by creating a proxy_summary data frame with the basic values for all the final proxies.

Then the site-by-site calibration interval statistics in the list “site_stats” are collected into a data frame and written to “Site_Cal_stats_Summary.csv.” These are the instrumental era proxy-by-proxy calibration statistics missing from the Pages 2k documentation. You can see the overall average R2 is only 0.16 and the average p is 0.165, these are poor correlations overall and there are few good correlations in the whole batch, the maximum R2 is only 0.68.

I should mention that a low calibration R² does not necessarily imply that a proxy is invalid. It may indicate that the proxy responds to seasonal, oceanographic, hydrologic, or regional variables not fully captured or represented by the nearest NOAA annual-mean temperature grid cell. Many of the proxies just represent one season which is not taken into account in Pages2k.r. Under some circumstances Pages 2k takes season into account but usually treats all proxies as if they were annual. The Pages 2K reconstruction also tries to take spatial distribution of the proxies into account, not something done by this program. Pages 2k also converts each series into a Z-score before merging them. I did try converting the proxies to Z-scores, as opposed to using the actual proxy measurements, and found no significant difference. With the noted exceptions, the Pages 2k NGEO supplement describes a procedure similar to that in Pages2k.R, except they used an infilled HadCRUT4 grid, whereas I use the NOAA GlobalTemp grid. Pages 2k does not provide proxy-by-proxy R2 values, but they do provide them region-by-region.

In the tropics, summer proxies are decent predictors of the annual mean, but in the higher latitudes they are not. In the extratropics and high latitudes annual mean temperature trends correlate to winter temperature trends (Pages 2K, 2017). One other point, the correlations used in Pages2k.r are all performed with ordinary least squares (the R function “lm(temp~proxy)”) which assumes no error in the proxy measurement. This is not always the case but here it was used for simplicity and as an exploratory tool.

Finally, the basic useful statistics for each proxy are written out to “Final_Proxy_List Summary.csv”, such as the proxy type, minimum and maximum values and the beginning and ending years. Counts of each proxy type (including those proxies with multiple measurements) are written out to “Final_Variable Count.csv.” Calibration interval correlation statistics are written to “Site_Cal_stats_Summary.csv.”

A programming tip. If you don’t understand a code block, cut-and-paste it into your favorite AI program and ask it to explain what the code does. Likewise, if you want to change what the program does or add a new code block, describe what you want to do to AI and ask it to write the new code for you. Be careful writing your coding request and you will probably get working code out of the AI, they can all write R code.

Published by Andy May

Petrophysicist, details available here: https://andymaypetrophysicist.com/about/

Leave a Reply