Load required packages.

library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.2
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(purrr)
## Warning: package 'purrr' was built under R version 3.5.2
library(readr)
library(tidyr)
library(ggplot2)
library(broom)
library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(truncnorm)

Composite indices of biomass burning

Define lakes for biomass burning analysis.

# Noatak River Basin (data from Higuera et al. 2011)
noatak <- c('LI','RA','PO','UC')  
# Brooks Range (data from Higuera et al. 2007,2009)
brooks <- c('CO','LC','RP','WK','XI') 
# Yukon Flats (data from Kelly et al. 2014)
yukon  <- c('CP','EP','GA','JA','LD','LT','LU','NR','PI','RE','RO','SL','WC','WI') 
# Copper River Basin (data from Barett et al. 2013)
copper <- c('CR','HD','MP1','SC') 

alaska <- c(noatak,brooks,yukon,copper)

Load functions to implement procedure: maximum likelihood estimation of the zero-inflated log-normal distributions in continuous 100-yr and 500-yr moving windows from 10,000 yr BP to -65 yr BP. Run for Alaska and all regions using specificed bandwidths (half the width of the moving window, within which the ZIL distribution will be estimated).

source(file.path('code','ZIL_Run.r'))

bandWidth   <- 50 
output100 <- map(list(alaska, noatak, brooks, yukon, copper),
                   zil_fn)

compCHAR_100 <- output100 %>% 
  `names<-` (c('Alaska','Noatak River Watershed','Kobuk Ridges and Valleys','Yukon Flats','Copper River Basin')) %>%
  bind_rows(.id = 'region')

# For reference, hard-coded for now
spans <- compCHAR_100 %>% 
  group_by(region) %>% 
  summarise(span = (1000/10) / n()) 

# Calculate ~500-year mean from loess regression
compositeCHAR <- compCHAR_100 %>% 
  group_by(region) %>% 
  do(augment(loess(composite.mean ~ yr.bp, span = 0.1, data = .), .)) %>% 
  rename(mean_100 = composite.mean,
         mean_500 = .fitted) %>% 
  select(region, yr.bp, mean_100, mean_500, CI.lower, CI.upper) %>% 
  do(augment(loess(CI.lower ~ yr.bp, span = 0.1, data = .), .)) %>%
  rename(CI_lower_500 = .fitted) %>%
  select(-.se.fit, -.resid) %>% 
  do(augment(loess(CI.upper ~ yr.bp, span = 0.1, data = .), .)) %>%
  rename(CI_upper_500 = .fitted) %>%
  select(-.se.fit, -.resid)

Quick plot of 100-year (thin line) and 500-year (thick line w/ CI) Alaska-wide and ecoregional composite indices of biomass burning.

ggplot(compositeCHAR) + 
    geom_line(aes(x= yr.bp, y = mean_100), color = 'grey40') +
    geom_ribbon(aes(x=yr.bp, ymin = CI_lower_500, ymax = CI_upper_500), alpha = 0.2) +
    geom_line(aes(x=yr.bp,y=mean_500), size = 1) +
    scale_x_reverse(breaks = seq(10000,-50,-2000)) +
    scale_y_continuous(breaks = c(1,2,3,4,5,6,7), labels =c(1,2,3,4,5,6,7)) +
    coord_cartesian(ylim = c(0.3,4)) +
    ylab("Biomass burning\n(standardized CHAR)") +
    xlab('Time (cal years BP)') +
    facet_wrap(~region) +
    theme_bw()

Re-create Figure 2, illustrating zero-inflated log-normal distribution of CHAR.

source(file.path('code','plotting_fig_2.r'))
zil_method_plot

Composite indices of fire frequency / fire synchrony

Percent sites burned per century

Redefine some regions (remove lakes with signal-to-noise index < 3)

## Update ecoregion parameters
#  MP1 in Copper River  and XI in Kobuk not suitable for peak analysis
copper <- c("CR","SC","HD")
brooks <- c("CO","LC","RP","WK")

# Update 'regions' object
regions  <- list('noatak'= noatak,'brooks'= brooks,'yukon'= yukon, 'copper'= copper) %>% 
  stack() %>% 
  rename(lake = values, region = ind)
alaska <- c(noatak,brooks,yukon,copper)

Load synchrony (percent of sites burned per century) analysis, as a function, stored in Synchrony_Run.r script.

source(file.path('code','synchrony_funs.r'))

Load the estimated timing of fire events from peak CHAR data.

fireYears_df <- fireyears_fn(alaska)

Run the regional version of the synchrony analysis. This result is presented in the supplemental informaton.

regional <- 1 #1 for yes, group by region, 0 for no, do all AK together
pctBurned_reg <- synch_fn(fireYears_df)

View those results (Figure S2):

Run the Alaska-wide version of the synchrony analysis.

regional <- 0 #1 for yes, group by region, 0 for no, do all AK together
pctBurned_ak <- synch_fn(fireYears_df)

pctBurned_ak <- pctBurned_ak %>% 
  augment(loess(boot.median ~ year, span = 0.1, data = .), .) %>% 
  rename(smooth_mean = .fitted) %>% 
  select(-.se.fit, -.resid) %>% 
  do(augment(loess(boot.lower ~ year, span = 0.1, data = .), .)) %>%
  rename(smooth_lower = .fitted) %>%
  select(-.se.fit, -.resid) %>% 
  do(augment(loess(boot.upper ~ year, span = 0.1, data = .), .)) %>%
  rename(smooth_upper = .fitted) %>%
  select(-.se.fit, -.resid)

View those results:

Fire return intervals

A 500-yr moving window mean estimate of FRI from peak-CHAR-estimated fire years was calcualted by PEH in MatLab.

FRIsmooth = read_csv(file.path('data','peakData','FRI_smooth.csv'))
FRIsmooth <- FRIsmooth %>% 
  filter(year_BP <= 8000)

Combine CHAR and FRI and calculate ratio of CHAR to FRI (i.e., fire severity)

ratio <- compositeCHAR %>% 
  filter(region == 'Alaska') %>% 
  right_join(FRIsmooth, by = c('yr.bp' = 'year_BP')) %>%
  mutate(ratio = mean_500/mFRI_1000,
         ratio.CIlower = CI_lower_500/CI_l,
         ratio.CIupper = CI_upper_500/CI_u) %>%
  select(yr.bp, ratio, ratio.CIlower, ratio.CIupper)

Climate data

Load climate data from Kaufman et al. (2017) and Wiles et al. (2014), and do some light munging.

# Wiles et al. 2014 data from the Gulf of Alaska ("GOA")
goa.df = read_csv(file.path('data','climateData','wiles_goa.csv'))
goa.mean <- mean(goa.df$temp)
goa.df$yearBP <- (1950 - goa.df$yearCE) 

goa.df$yearBins <- cut(goa.df$yearBP,include.lowest = T,right = F,
                       breaks = seq(min(goa.df$yearBP),max(goa.df$yearCE),10), 
                       labels = seq(min(goa.df$yearBP),2000,10))
goa.df$zscore = (goa.df$temp - goa.mean) 

goa.binned <- goa.df %>%
  group_by(yearBins) %>%
  summarise(bin.temp = mean(temp)) %>%
  mutate(zscore = (bin.temp - goa.mean)) %>%
  mutate(sign = ifelse(zscore >= 0,'positive','negative'))

# Kaufman et al. pan-Alaska dataset 
clim.df = read_csv(file.path('data','climateData','kaufman_paleoclimate_synthesis.csv'))
clim.df[["tempsign"]] <- ifelse(clim.df[["all.temp"]] >= 0, "positive", "negative")
clim.df[["moistsign"]] <- ifelse(clim.df[["all.moist"]] >= 0, "positive", "negative")

Re-create Figure 5

source(file.path('code','plotting_fig_suppver5.r'))
fig_5

Correlation analysis

Read in climate data from Kaufman et al. XXXX, a composite Holocene paleoclimate record for Beringia

kaufman_df <- read_csv(file.path('data','climateData','kaufman_paleoclimate_synthesis.csv'))

Last 10,000 years

Join all of the charcoal-based metrics together with the climate record, binned into 500-year intervals that correspod with Kaufman record.

corr_df_10kya <- compositeCHAR %>% 
  filter(region == 'Alaska') %>% 
  right_join(FRIsmooth, by = c('yr.bp' = 'year_BP')) %>%
  right_join(ratio, by = 'yr.bp') %>%
  full_join(pctBurned_ak, by = c('yr.bp' = 'year')) %>%
  ungroup() %>% 
  dplyr::select(-region.x) %>% 
  mutate(bins500 = as.numeric(as.character(cut(.$yr.bp,include.lowest = T, right = T,
                                            breaks = seq(0,10000,500), 
                                            labels = seq(250,10000,500)))),
         bins500 = ifelse(yr.bp < 0, 250, bins500)) %>%
  group_by(bins500) %>%
  summarise_if(is.numeric, mean, na.rm = T) %>%
  select(-yr.bp) %>% 
  rename(yr.bp = bins500) %>%
  full_join(kaufman_df, by = 'yr.bp') %>%
  filter(yr.bp <= 8000) %>%
  dplyr::select(CHAR_500 = mean_500, 
         CHAR_100 = mean_100, 
         Pct_Burned = boot.median, 
         mFRI_1000, 
         ratio,
         all_temp = all.temp,
         all_moist = all.moist,
         midgeA = midgeA.temp,
         midgeB = midgeB.temp,
         pollenAnnual = pollenAnn.temp,
         pollenSummer = pollenSumm.temp)

Calculate Pearson pairwise correlations.

Hmisc::rcorr(as.matrix(corr_df_10kya), type="pearson") # type can be pearson or spearman
##              CHAR_500 CHAR_100 Pct_Burned mFRI_1000 ratio all_temp
## CHAR_500         1.00     1.00       0.46     -0.60  0.96    -0.27
## CHAR_100         1.00     1.00       0.48     -0.61  0.95    -0.26
## Pct_Burned       0.46     0.48       1.00     -0.92  0.68    -0.77
## mFRI_1000       -0.60    -0.61      -0.92      1.00 -0.80     0.79
## ratio            0.96     0.95       0.68     -0.80  1.00    -0.52
## all_temp        -0.27    -0.26      -0.77      0.79 -0.52     1.00
## all_moist        0.21     0.23       0.29     -0.30  0.24    -0.02
## midgeA          -0.09    -0.09      -0.38      0.36 -0.24     0.69
## midgeB          -0.34    -0.35      -0.73      0.65 -0.50     0.67
## pollenAnnual    -0.14    -0.11      -0.35      0.36 -0.27     0.60
## pollenSummer    -0.07    -0.06      -0.14      0.18 -0.13     0.19
##              all_moist midgeA midgeB pollenAnnual pollenSummer
## CHAR_500          0.21  -0.09  -0.34        -0.14        -0.07
## CHAR_100          0.23  -0.09  -0.35        -0.11        -0.06
## Pct_Burned        0.29  -0.38  -0.73        -0.35        -0.14
## mFRI_1000        -0.30   0.36   0.65         0.36         0.18
## ratio             0.24  -0.24  -0.50        -0.27        -0.13
## all_temp         -0.02   0.69   0.67         0.60         0.19
## all_moist         1.00   0.10  -0.32         0.15        -0.31
## midgeA            0.10   1.00   0.74         0.52         0.35
## midgeB           -0.32   0.74   1.00         0.30         0.22
## pollenAnnual      0.15   0.52   0.30         1.00         0.67
## pollenSummer     -0.31   0.35   0.22         0.67         1.00
## 
## n= 16 
## 
## 
## P
##              CHAR_500 CHAR_100 Pct_Burned mFRI_1000 ratio  all_temp
## CHAR_500              0.0000   0.0713     0.0132    0.0000 0.3053  
## CHAR_100     0.0000            0.0624     0.0120    0.0000 0.3290  
## Pct_Burned   0.0713   0.0624              0.0000    0.0041 0.0005  
## mFRI_1000    0.0132   0.0120   0.0000               0.0002 0.0003  
## ratio        0.0000   0.0000   0.0041     0.0002           0.0410  
## all_temp     0.3053   0.3290   0.0005     0.0003    0.0410         
## all_moist    0.4337   0.4018   0.2745     0.2579    0.3776 0.9391  
## midgeA       0.7475   0.7442   0.1457     0.1693    0.3774 0.0032  
## midgeB       0.1964   0.1817   0.0014     0.0067    0.0497 0.0043  
## pollenAnnual 0.6119   0.6922   0.1776     0.1735    0.3074 0.0144  
## pollenSummer 0.7976   0.8290   0.6003     0.5038    0.6333 0.4850  
##              all_moist midgeA midgeB pollenAnnual pollenSummer
## CHAR_500     0.4337    0.7475 0.1964 0.6119       0.7976      
## CHAR_100     0.4018    0.7442 0.1817 0.6922       0.8290      
## Pct_Burned   0.2745    0.1457 0.0014 0.1776       0.6003      
## mFRI_1000    0.2579    0.1693 0.0067 0.1735       0.5038      
## ratio        0.3776    0.3774 0.0497 0.3074       0.6333      
## all_temp     0.9391    0.0032 0.0043 0.0144       0.4850      
## all_moist              0.7163 0.2256 0.5905       0.2403      
## midgeA       0.7163           0.0011 0.0369       0.1811      
## midgeB       0.2256    0.0011        0.2568       0.4037      
## pollenAnnual 0.5905    0.0369 0.2568              0.0046      
## pollenSummer 0.2403    0.1811 0.4037 0.0046

Last 1200 years

wiles_df <- read_csv(file.path('data','climateData','wiles_goa.csv')) %>% 
  mutate(yr.bp = 1950-yearCE)

corr_df_1kya <- wiles_df %>%
  mutate(bins10 = rev(cut(.$yr.bp,include.lowest = T, right = T,
                          breaks = seq(1150,-60,-10),
                          labels = seq(1150,-50,-10)))) %>%
  mutate(bins10 = as.numeric(as.character(bins10))) %>% 
  group_by(bins10) %>%
  summarise_all(.funs = mean) %>%
  select(-yr.bp, -yearCE) %>%
  rename(yr.bp = bins10) %>%
  full_join(.,filter(compositeCHAR, region == 'Alaska')) %>%
  full_join(pctBurned_ak, by = c('yr.bp' = 'year')) %>%
  full_join(FRIsmooth, by = c('yr.bp' = 'year_BP')) %>%
  full_join(ratio) %>% 
  select(CHAR_500 = mean_500, 
         CHAR_100 = mean_100, 
         Pct_Burned = boot.median, 
         mFRI_1000, 
         ratio,
         GOAtemp = temp)

Calculate Pearson pairwise correlations.

Hmisc::rcorr(as.matrix(corr_df_1kya), type="pearson") # type can be pearson or spearman
##            CHAR_500 CHAR_100 Pct_Burned mFRI_1000 ratio GOAtemp
## CHAR_500       1.00     0.95       0.32     -0.47  0.97    0.43
## CHAR_100       0.95     1.00       0.35     -0.43  0.90    0.41
## Pct_Burned     0.32     0.35       1.00     -0.65  0.50    0.37
## mFRI_1000     -0.47    -0.43      -0.65      1.00 -0.66   -0.41
## ratio          0.97     0.90       0.50     -0.66  1.00    0.43
## GOAtemp        0.43     0.41       0.37     -0.41  0.43    1.00
## 
## n
##            CHAR_500 CHAR_100 Pct_Burned mFRI_1000 ratio GOAtemp
## CHAR_500       1007     1007        803       162   162     121
## CHAR_100       1007     1007        803       162   162     121
## Pct_Burned      803      803        803       161   161     118
## mFRI_1000       162      162        161       162   162      24
## ratio           162      162        161       162   162      24
## GOAtemp         121      121        118        24    24     121
## 
## P
##            CHAR_500 CHAR_100 Pct_Burned mFRI_1000 ratio  GOAtemp
## CHAR_500            0.0000   0.0000     0.0000    0.0000 0.0000 
## CHAR_100   0.0000            0.0000     0.0000    0.0000 0.0000 
## Pct_Burned 0.0000   0.0000              0.0000    0.0000 0.0000 
## mFRI_1000  0.0000   0.0000   0.0000               0.0000 0.0453 
## ratio      0.0000   0.0000   0.0000     0.0000           0.0375 
## GOAtemp    0.0000   0.0000   0.0000     0.0453    0.0375