Forest manipulation experiment reveals divergent controls on the sources and age of lateral DOC and CO₂ export

Authors
Affiliation

Audrey Campeau

Université de Montréal & SLU

A. Zannella and M. Wallin

Published

August 22, 2025

1 Introduction

1.1 General Context

  • Lateral C export is a significant fraction of watershed C balance.

  • Forested catchments contain a large OM storage that could sustain DOC export for decades (Ledesma et al. 2013)

  • LCE and NEE are connected over long timescale, by hydrology (Öquist et al. 2014)

  • Isotopic values of C (stable and radiogenic) can inform us on the sources and age of C.

1.2 Research Question

  • What are the controls over the sources and age of lateral CO2 and DOC export in forested catchments?

  • Can a forest manipulation experiment (forest clearcut and ditch cleaning) provide new insight to test ongoing hypothesis on the controls of LCE in forested catchments?

1.3 Hypothesis

  • The CO2 source and age is more closely linked to the forest C sink (A. Campeau et al. 2019), so clearcutting the forest should have an impact on C sources and age

  • The DOC source and age is linked to discharge (Audrey Campeau et al. 2017) or water table position (A. Campeau et al. 2019), so changes in watershed hydrology, caused by clear-cutting and draining, should change the source and age of DOC.

1.4 Main Conclusion

DOC is controled more by hydrological processes, which determines what material is being mobilised, while CO2 is controled more by biological processes, which fuels CO2 in the watershed. Both are therefore controled by different processes, but will likely respond to changes in climate, albeit via different drivers.

2 Methodology

2.1 Study Site and Treatment:

  • Six headwater catchments are included in this study:

    • 2 pristine sites (C1 and C2)

    • 4 treated watersheds (DC1 to DC4).

  • The DC sites received different treatments:

    • Forest in all four sites was clearcut - around July 2020.

    • Two sites, DC1 and DC3, were also ditch cleaned - in September 2021.

    • The treatments are named as follow (pristine, clear-cut and ditch cleaning)

2.2 Field measurements:

  • All four sites are monitored for flow and water chemistry on a near continuous basis.

  • Radiocarbon and stable C isotope measurements were collected at those six sites simultaniously and throughout various treatment stages.

  • 14C measurements

    • Start 2020-03-12

    • End 2022-10-25

2.3 Statistical analysis

  • Dunn’s test for non-parametric grouping

  • Linear regression models, combined with ANCOVA to identify significant differences in intercept or slope.

  • Linear mixed-effect models would be possible if we wanted to build a predictive model, but this is not the aim here.

3 Results

3.1 Hydrograph, C concentrations and ¹⁴C-content across sites over time

👈 see code here 
# Define operation periods for treatments::::::::::::::::::::::::::::::::::::::::::
  
clearcut_start <- as.Date("2020-07-01")
clearcut_end <- as.Date("2020-08-25")
ditch_cleaning_start <- as.Date("2021-09-01")
ditch_cleaning_end <- as.Date("2021-09-30")

 
# Make the precipitation graph :::::::::::::::::::::::::::::::::::::::::::::::::::::


precipitation= ggplot(DC_Q %>%
         filter(Site_id %in% c("DC1")),
       aes(y = P, x = as.Date(Date), color = Site_id)) +
  
  # Add gray background bars for operations
  annotate("rect", 
           xmin = clearcut_start, xmax = clearcut_end,
           ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3) +
  
  annotate("rect", 
           xmin = ditch_cleaning_start, xmax = ditch_cleaning_end,
           ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3) +
  
  geom_bar(stat = 'identity', width=0.2) +
  labs( y= "Precipitation (mm/day)") +
  scale_y_reverse()+
  scale_color_manual(values=site_colors_6)+ #"
  scale_x_date(limits=c(as.Date("2020-01-01"), as.Date("2022-11-01")),
               date_labels = "%b-%Y", date_breaks = "6 months") +
  
  
  labs(y = "P (mm/d)",  x = "Date") +
  theme(legend.position = "none", axis.title.x = element_blank())



# Make Hydrograph :::::::::::::::::::::::::::::::::::::::::::::::::::::

hydrograph= ggplot(DC_Q %>%
         filter(Site_id %in% c("DC2", "DC3")),
       aes(y = q_md*1000, x = as.Date(Date), color = Site_id)) +
  
  # Add gray background bars for operations
  annotate("rect", 
           xmin = clearcut_start, xmax = clearcut_end,
           ymin = -Inf, ymax = Inf,
           fill = "gray70", alpha = 0.3) +
  
  annotate("rect", 
           xmin = ditch_cleaning_start, xmax = ditch_cleaning_end,
           ymin = -Inf, ymax = Inf,
           fill = "gray70", alpha = 0.3) +
  
  # Add the discharge lines
  geom_line() +
  scale_color_manual(values=site_colors_6)+ #"
   
  # Add text labels for operations
 annotate("text", x = clearcut_start + (clearcut_end - clearcut_start)/2, y = Inf, 
           label = "Forest clear cut operations\n(July 2020)",
           vjust = 1.2, hjust = 0.5, size = 3.5, color = "black") +
  
  annotate("text", x = ditch_cleaning_start + (ditch_cleaning_end -    ditch_cleaning_start)/2, y = Inf, 
           label = "Ditch cleaning operations\n(September 2021)",
           vjust = 1.2, hjust = 0.5, size = 3.5, color = "black") +
  
  scale_x_date(limits=c(as.Date("2020-01-01"), as.Date("2022-11-01")),
                 date_labels = "%b-%Y", date_breaks = "6 months") +
  labs(y = "q (mm/d)",   x = "Date",  color = "Stream") +
  theme(legend.position = "none", axis.title.x = element_blank())



ggarrange(precipitation, hydrograph,
          ncol=1, align="v")

Precipitation and discharge timeseries
👈 see code here 
# Timeseries 14C per site, CO2 and DOC seperately
plot_timeseries = function(data, y_var,labs_y) {
  ggplot(data, aes(x=Date, y=!!sym(y_var), group=Site_id, color=Site_id)) +
    geom_line(data=data[!is.na(data[[y_var]]), ]) +
    geom_point(size=3) +
  
  # Add gray background bars for operations
  annotate("rect", 
           xmin = clearcut_start, xmax = clearcut_end,
           ymin = -Inf, ymax = Inf,
           fill = "gray70", alpha = 0.3) +
  
  annotate("rect", 
           xmin = ditch_cleaning_start, xmax = ditch_cleaning_end,
           ymin = -Inf, ymax = Inf,
           fill = "gray70", alpha = 0.3) +
    
    scale_x_date(limits=c(as.Date("2020-01-01"), as.Date("2022-11-01")),
                 date_labels = "%b-%Y", date_breaks = "6 months") +
    scale_color_manual(values=c(site_colors_6))+
    labs(y=labs_y)+
    
    theme( # set margins to align combined plots 
      plot.margin = margin(5.5, 5.5, 5.5, 5.5, "pt"),
      axis.text.x = element_text(size = 10),
      axis.title.x = element_text(size = 11),
      panel.grid.major = element_line(size = 0.5),
      panel.grid.minor = element_line(size = 0.25)
    )

}



#Concentration timeseries


timeserie_DOC = ggExtra::ggMarginal(
  plot_timeseries(DC_Q, "DOC_mgL", "DOC (mg/L)")+
    # Add text labels for operations
 annotate("text", x = clearcut_start + (clearcut_end - clearcut_start)/2, y = Inf, 
           label = "Forest clear cut operations\n(July 2020)",
           vjust = 1.2, hjust = 0.5, size = 3.5, color = "black") +
  
  annotate("text", x = ditch_cleaning_start + (ditch_cleaning_end -    ditch_cleaning_start)/2, y = Inf, 
           label = "Ditch cleaning operations\n(September 2021)",
           vjust = 1.2, hjust = 0.5, size = 3.5, color = "black") +
  theme(legend.position = "top", axis.title.x = element_blank()),
  margins=c("y"), type = "boxplot", groupColour = TRUE, groupFill = TRUE )


timeserie_14CDOC = ggExtra::ggMarginal(
  plot_timeseries(DC_Q, "DOC_14C_Modern", "14C-DOC (%)")+
  scale_y_continuous(limits=c(95,115))+
  theme(legend.position = "none", axis.title.x = element_blank()),
  margins=c("y"), type = "boxplot", groupColour = TRUE, groupFill = TRUE )



timeserie_CO2 = ggExtra::ggMarginal(
  plot_timeseries(DC_Q, "CO2_mgL", "CO2 (mg/L)")+
     theme(legend.position = "none", axis.title.x = element_blank()),
  margins=c("y"), type = "boxplot", groupColour = TRUE, groupFill = TRUE )

timeserie_14CCO2 = ggExtra::ggMarginal(
  plot_timeseries(DC_Q, "CO2_14C_Modern", "14C-CO2 (%)")+
  scale_y_continuous(limits=c(95,115))+
   theme(legend.position = "none", axis.title.x = element_blank()),
  margins=c("y"), type = "boxplot", groupColour = TRUE, groupFill = TRUE )


# Combine four plts
library(gridExtra)



final_plot <- grid.arrange(
  timeserie_DOC,
  timeserie_14CDOC,
  timeserie_CO2,
  timeserie_14CCO2,
  ncol = 1,
  heights = c( 0.75, 0.5, 0.5, 0.55)
)

C concentration and radiocarbon content timeseries coloured by sites

3.1.0.1 Dunn’s test | ¹⁴C and [C] ~ Site

👈 see code here 
# Calculate medians for each variable by Site_id
combined_medians <- DC_Q %>%
  group_by(Site_id) %>%
  summarise(
    "[DOC]" = round(median(DOC_mgL, na.rm = TRUE), 2),
    "14C-DOC" = round(median(DOC_14C_Modern, na.rm = TRUE),2),
    "[CO2]" = round(median(CO2_mgL, na.rm = TRUE), 2),
    "14C-CO2" = round(median(CO2_14C_Modern, na.rm = TRUE),2),
    .groups = 'drop'
  ) %>%
  rename(Site = Site_id)

knitr::kable(combined_medians, 
             caption = "Median by site for DOC and CO2 concentration and 14C-content")
Median by site for DOC and CO2 concentration and 14C-content
Site [DOC] 14C-DOC [CO2] 14C-CO2
DC2 39.14 106.82 2.17 104.01
DC4 44.10 106.70 4.43 101.54
DC1 21.40 106.74 1.66 104.47
DC3 22.05 105.86 1.26 102.99
👈 see code here 
library(dunn.test)
library(multcompView)
library(PMCMR)
# Perform the test and get multcomp letters
DOCSites_letters = multcompLetters(get.pvalues(PMCMRplus::kwAllPairsDunnTest(
          DC_Q$DOC_mgL ~ DC_Q$Site_id, p.adjust="bonf")),
          threshold=0.05)$Letters

DOC_14C_Sites_letters = multcompLetters(get.pvalues(PMCMRplus::kwAllPairsDunnTest(
          DC_Q$DOC_14C_Modern ~ DC_Q$Site_id, p.adjust="bonf")),
          threshold=0.05)$Letters
        
CO2Sites_letters = multcompLetters(get.pvalues(PMCMRplus::kwAllPairsDunnTest(
          DC_Q$CO2_mgL ~ DC_Q$Site_id, p.adjust="bonf")),
          threshold=0.05)$Letters

CO2_C14_Sites_letters = multcompLetters(get.pvalues(PMCMRplus::kwAllPairsDunnTest(
          DC_Q$CO2_14C_Modern ~ DC_Q$Site_id, p.adjust="bonf")),
          threshold=0.05)$Letters

        
# Combine into a data frame
combined_dunns_results <- data.frame(
  Site = names(DOCSites_letters),
  "[DOC]" = as.character(DOCSites_letters),
  "14C-DOC" = as.character(DOC_14C_Sites_letters),
  "[CO2]" = as.character(CO2Sites_letters[names(DOCSites_letters)]),
  "14C-CO2" = as.character(CO2_C14_Sites_letters[names(DOCSites_letters)]),
  check.names = FALSE
)

knitr::kable(combined_dunns_results, 
             caption = "Grouping letters for dunn's test result")
Grouping letters for dunn’s test result
Site [DOC] 14C-DOC [CO2] 14C-CO2
DC2 a a a ab
DC4 a a b a
DC1 b a ac b
DC3 b a c ab
Interpretation

Trend:

  • No obvious trends in C concentration or 14C content over time

  • There was a clear peak in DOC and CO2 concentration at DC4 during clearcut opperations.

Differences between sites:

  • CO2 and DOC concentration at DC4 is consistently higher than the other sites, followed with DC2, and DC1 and DC3.

  • 14C-DOC is not different across sites

  • 14C-CO2 is significantly lower at DC4, followed with DC3 and DC2, DC1 which is significantly higher

Q data

The database contains DC3 and DC2 flow data, one ditch cleaned while the other only clearcut. Other timeseries are incomplete.

3.2 Changes in C concentration or ¹⁴C-content across Treatments

Is there a significant change in the median 14C content of CO2 and DOC between sites or treatment, based on their distribution ?

👈 see code here 
# Make a scatterplot to distinguish the grouping across sites
density_14Cscatter_sites = 
  ggExtra::ggMarginal(
  
            ggplot(data=DC_Q,
                 aes(x=DOC_mgL, y=CO2_mgL_filled, 
                    color=Treatment))+
            stat_density_2d(aes(fill = Treatment), geom = "polygon", linewidth=0,
                             alpha = 0.2)+ 
            geom_point(size=1)+
            geom_abline(intercept=0,slope=1, linetype="dotted")+
            geom_smooth(method="glm", se=F, aes(color=Treatment), show.legend = T)+
  
            stat_regline_equation(label.y.npc = "top", label.x.npc = "left",
            aes(label =  paste(..eq.label.., ..adj.rr.label.., 
                               sep = "~~~~"), color = Treatment), 
                                show.legend = F, size=3)+
    
            scale_color_manual(values=c(treatments_colors))+
            scale_fill_manual(values=c(treatments_colors))+
            scale_y_log10()+
              scale_x_log10()+
            
            labs(x=bquote("DOC (mg L"^-1*")"), 
                   y=bquote("CO"[2]* " (mg L"^-1*")"))+
              theme(legend.position = "none"), 
            type = "boxplot", groupColour = TRUE, groupFill = TRUE ) 


# Make a scatterplot to distinguish the grouping across treatment
density_14Cscatter_treatment = 
ggExtra::ggMarginal(
  
            ggplot(data=DC_Q,aes(x=DOC_14C_Modern, y=CO2_14C_Modern, 
                    color=Treatment))+
               stat_density_2d(aes(fill = Treatment), geom = "polygon", linewidth=0,
                             alpha = 0.2)+ 
            geom_point(size=2)+
            geom_abline(intercept=0,slope=1, linetype="dotted")+
            scale_color_manual(values=c(treatments_colors))+
            scale_fill_manual(values=c(treatments_colors))+
            scale_y_continuous(limit=c(95,115))+
            labs(x=bquote(Delta^14*"C-DOC (% modern)"), 
                   y=bquote(Delta^14*"C-CO"[2]*" (% modern)"))+
              theme(legend.position = "bottom"), 
            type = "boxplot", groupColour = TRUE, groupFill = TRUE )
            


ggarrange(density_14Cscatter_sites,
density_14Cscatter_treatment, nrow=2, align="hv")

relationship between C concentration and radiocarbon content of CO2 and DOC coloured by treatment

3.2.0.1 ANCOVA | [CO2] ~ [DOC] / Treatment

👈 see code here 
#Run ANCOVA test - Difference in hydro. response between treatments


knitr::kable(DC_Q  %>%  anova_test(DOC_mgL ~ CO2_mgL_filled*Treatment), 
             caption = "Ancova test results: DOC ~ CO2 x Treatment")
Ancova test results: DOC ~ CO2 x Treatment
Effect DFn DFd F p p<.05 ges
CO2_mgL_filled 1 298 526.166 0e+00 * 0.638
Treatment 2 298 28.855 0e+00 * 0.162
CO2_mgL_filled:Treatment 2 298 16.675 1e-07 * 0.101

3.2.0.2 Dunn’s test | ¹⁴C et [C]~ Treatment

👈 see code here 
# Calculate medians for each variable by Site_id
combined_medians <- DC_Q %>%
  group_by(Treatment) %>%
  summarise(
    "14C-DOC" = round(median(DOC_14C_Modern, na.rm = TRUE),2),
     "[DOC]" = round(median(DOC_mgL, na.rm = TRUE), 2),
    "14C-CO2" = round(median(CO2_14C_Modern, na.rm = TRUE),2),
    "[CO2]" = round(median(CO2_mgL, na.rm = TRUE), 2),
    .groups = 'drop'
  ) %>%
  rename(Site = Treatment)

knitr::kable(combined_medians, 
             caption = "Median by treatment for DOC and CO2 concentration and 14C-content")
Median by treatment for DOC and CO2 concentration and 14C-content
Site 14C-DOC [DOC] 14C-CO2 [CO2]
Pristine 105.53 26.48 101.06 1.95
Clearcut 106.90 39.28 103.49 2.15
Ditch cleaning 106.49 21.15 104.30 1.42
👈 see code here 
library(dunn.test)
library(multcompView)
library(PMCMR)

# Perform the test and get multcomp letters
DOCTreatment_letters = multcompLetters(get.pvalues(PMCMRplus::kwAllPairsDunnTest(
          DC_Q$DOC_mgL ~ DC_Q$Treatment, p.adjust="bonf")),
          threshold=0.05)$Letters
        
DOC_C14_Treatment_letters = multcompLetters(get.pvalues(PMCMRplus::kwAllPairsDunnTest(
          DC_Q$DOC_14C_Modern ~ DC_Q$Treatment, p.adjust="bonf")),
          threshold=0.05)$Letters
        
CO2Treatment_letters = multcompLetters(get.pvalues(PMCMRplus::kwAllPairsDunnTest(
          DC_Q$CO2_mgL ~ DC_Q$Treatment, p.adjust="bonf")),
          threshold=0.05)$Letters

CO2_C14_Treatment_letters = multcompLetters(get.pvalues(PMCMRplus::kwAllPairsDunnTest(
          DC_Q$CO2_14C_Modern ~ DC_Q$Treatment, p.adjust="bonf")),
          threshold=0.05)$Letters

        
# Combine into a data frame
combined_dunns_results <- data.frame(
  Site = names(DOC_C14_Treatment_letters),
  "14C-DOC" = as.character(DOC_C14_Treatment_letters),
  "[DOC]" = as.character(DOCTreatment_letters),
  "14C-CO2" = as.character(CO2_C14_Treatment_letters[names(DOC_C14_Treatment_letters)]),
  "[CO2]" = as.character(CO2Treatment_letters[names(DOC_C14_Treatment_letters)]),
  check.names = FALSE
)


knitr::kable(combined_dunns_results, 
             caption = "Grouping letters for dunns test C concentration and 14C content by treatment")
Grouping letters for dunns test C concentration and 14C content by treatment
Site 14C-DOC [DOC] 14C-CO2 [CO2]
Pristine a a a a
Clearcut b b a a
Ditch cleaning ab a a a
Interpretation

Relationship between CO2 and DOC concentration:

  • The relationship between CO2 and DOC concentration is positive for pristine condition, but the slope becomes less significant with clearcut (closer to 0) and not significant (slope =0) with ditch cleaning. Hence the DOC concentration continues to vary across a wide range with treatments, but the CO2 becomes more stable and independant of DOC concentration. This could be due to changes in their mobilization pathways or simply because clearcuting “removes” a CO2 source.

Treatment effect on C concentrations:

  • The DOC concentrations are significantly higher following clearcut treatment compared with ditchcleaned and pristine conditions (short term effect of ditch cleaning can compensate?)

  • The CO2 concentration do not differ significantly across treatment

Treatment effect on 14 content

  • The 14C-DOC is significantly lower in the pristine (group a) compared with clearcut and ditch cleaning sites.

  • The 14C-CO2 doesn’t differ significantly across treatments

3.3 Relationships - controls on C sources, age and concentrations

3.3.1 Hydrological control over C concentrations

Is the radiocarbon age or concentration of DOC and CO2 controled by runoff, and does this relationship changes following treatment?

👈 see code here 
# Hydrological control over  C concentration by Site
hydro_resp_C_site=
  ggplot(DC_Q %>% pivot_longer(
                     cols = c(DOC_mgL, CO2_mgL),
                     names_to = "Carbon_species",
                     values_to = "C_mgL"
                   ),
       aes(y=C_mgL, x=q_md_filled*1000, color=Site_id))+ 
  geom_point(size=1, show.legend = F)+
  facet_wrap(~Carbon_species)+
  scale_x_log10()+
  scale_y_log10()+
  geom_smooth(method="lm", se=F, aes(color=Site_id), show.legend = T)+
  scale_color_manual(values=site_colors_6)+
  stat_regline_equation(label.y.npc = "top", label.x.npc = "left",
  aes(label =  paste(..eq.label.., ..adj.rr.label.., 
                     sep = "~~~~"), color = Site_id), 
                      show.legend = F, size=3)+
  labs(x="q (mm/d)", 
       y=bquote("C (mg/L)"))+
  theme(legend.position = "top", base_size=8)


# Hydrological control over C concentration by treatment
hydro_resp_C_treatment=
  ggplot(DC_Q %>% pivot_longer(
                     cols = c(DOC_mgL, CO2_mgL),
                     names_to = "Carbon_species",
                     values_to = "C_mgL"
                   ),
       aes(y=C_mgL, x=q_md_filled*1000, color=Treatment))+ 
  geom_point(size=1, show.legend = F)+
  facet_wrap(~Carbon_species)+
  scale_x_log10()+
  scale_y_log10()+
  geom_smooth(method="lm", se=F, aes(color=Treatment), show.legend = T)+
  scale_color_manual(values=treatments_colors)+
  stat_regline_equation(label.y.npc = "top", label.x.npc = "left",
  aes(label =  paste(..eq.label.., ..adj.rr.label.., 
                     sep = "~~~~"), color = Treatment), 
                      show.legend = F, size=3)+
  labs(x="q (mm/d)", 
       y=bquote("C (mg/L)"))+
  theme(legend.position = "top", base_size=8)





ggarrange(hydro_resp_C_site, 
          hydro_resp_C_treatment, 
          nrow=2, common.legend = F)

relationship between C concetration and discharge coloured by site and treatment

3.3.1.1 ANCOVA test

Is there a significant difference in the hydrological response of DOC concentrations between Treatments or Sites?

👈 see code here 
#Run ANCOVA test - Difference in hydro. response between treatments

knitr::kable(DC_Q  %>%  anova_test(DOC_mgL ~ q_md_filled*Treatment), 
             caption = "ANCOVA test results - DOC * q / treatment")
ANCOVA test results - DOC * q / treatment
Effect DFn DFd F p p<.05 ges
q_md_filled 1 285 3.461 0.064 0.012
Treatment 2 285 18.080 0.000 * 0.113
q_md_filled:Treatment 2 285 1.277 0.280 0.009
👈 see code here 
knitr::kable(DC_Q  %>%  anova_test(DOC_mgL ~ q_md_filled*Site_id), 
             caption = "ANCOVA test results - DOC * q / Site")
ANCOVA test results - DOC * q / Site
Effect DFn DFd F p p<.05 ges
q_md_filled 1 283 0.002 9.66e-01 6.60e-06
Site_id 3 283 39.524 0.00e+00 * 2.95e-01
q_md_filled:Site_id 3 283 8.764 1.41e-05 * 8.50e-02

Is there a significant difference in the hydrological response of CO2 concentrations between Treatments or Sites?

👈 see code here 
#Run ANCOVA test - Difference in hydro. response between treatments
knitr::kable(DC_Q  %>%  anova_test(CO2_mgL ~ q_md_filled*Site_id), 
             caption = "ANCOVA test results - CO2 * q / treatment")
ANCOVA test results - CO2 * q / treatment
Effect DFn DFd F p p<.05 ges
q_md_filled 1 95 2.113 0.149 0.022
Site_id 3 95 31.327 0.000 * 0.497
q_md_filled:Site_id 3 95 3.818 0.012 * 0.108
👈 see code here 
knitr::kable(DC_Q  %>%  anova_test(CO2_mgL ~ q_md_filled*Site_id), 
             caption = "ANCOVA test results - CO2 * q / treatment")
ANCOVA test results - CO2 * q / treatment
Effect DFn DFd F p p<.05 ges
q_md_filled 1 95 2.113 0.149 0.022
Site_id 3 95 31.327 0.000 * 0.497
q_md_filled:Site_id 3 95 3.818 0.012 * 0.108
Interpretation

Hydrological response:

  • DOC and CO2 concentrations are not controled by discharge, except for DC4 where we believe there is a wetland source that could sustain flow in lowflow conditions.

  • ANCOVA indicates that Q doesn’t have a significant effect on DOC concentration, but treatment and site do (intercept difference). Significant interaction between Q and Site (slope difference)

  • ANCOVA indicates that Q doesnt have a significant effect on CO2 concentration, nor does treatment. Site effect causes a signficant effect on intercept and slope (interaction)

3.3.2 Hydrological controls over 14C-content

👈 see code here 
#Scatter coloured by Treatment
ggplot( DC_Q %>% pivot_longer(
                     cols = c(CO2_14C_Modern, DOC_14C_Modern),
                     names_to = "Carbon_specie",
                     values_to = "C14_value"
                   ),
       aes(y=C14_value, x=q_md_filled*1000, fill=Treatment, color=Treatment))+
  geom_point(size=3)+
  geom_smooth(method="glm", se=F, aes(color=Treatment), show.legend = F)+
  scale_fill_manual(values=treatments_colors)+ 
  scale_color_manual(values=treatments_colors)+ 

  stat_regline_equation(
  label.y.npc = "top", label.x.npc = "left",
  aes(label =  paste(..eq.label.., ..adj.rr.label.., 
                     sep = "~~~~"), color = Treatment), 
  show.legend = F, size=4)+
  labs(x="q (m/d)", y=bquote(Delta^14*"C (% modern)"), shape="Watershed ID")+
  facet_wrap(~Carbon_specie, scale="fixed", nrow=1)+
  theme(legend.position = "top")

relationship between radiocarbon content of CO2 and DOC and discharge, coloured by treatment

3.3.2.1 ANCOVA test

Is there a significant difference in the hydrological response of 14C-DOC between Treatments?

Ancova test results: 14C-DOC ~ q x Treatment
Effect DFn DFd F p p<.05 ges
q_md_filled 1 40 22.622 2.57e-05 * 0.361
Treatment 2 40 13.505 3.30e-05 * 0.403
q_md_filled:Treatment 2 40 2.234 1.20e-01 0.100

Is there a significant difference in the hydrological response of 14C-DOC between Sites?

ANCOVA test results, 14C-DOC * q / Site
Effect DFn DFd F p p<.05 ges
q_md_filled 1 38 11.645 0.002 * 0.235
Site_id 3 38 2.254 0.098 0.151
q_md_filled:Site_id 3 38 1.076 0.371 0.078

Is there a significant difference in the hydrological response of 14C-CO2 between Treatments?

Ancova test results: 14C-CO2 ~ q x Treatment
Effect DFn DFd F p p<.05 ges
q_md_filled 1 37 0.119 0.732 0.003
Treatment 2 37 0.677 0.514 0.035
q_md_filled:Treatment 1 37 0.282 0.598 0.008
Interpretation

  • While q doesn’t have a significant effect on DOC concentration, it does have a signficant effect on 14C-DOC. The 14C-content generally increases with Q, mobilizing more post-bomb C stored in soils.

  • Treatment has a significant effect in the ANCOVA, indicating that the intercept shift in the 14C-DOC~q relationship, from 100%modern in pristine sites to 110%modern in clearcut+ditchcleaned sites, is significant. The effect of Site is not significant.

  • Q doesnt have en effect but not on 14C-CO2

3.4 Biological controls over 14C-CO2 - Keeling plots

👈 see code here 
DC_Q$CO2_mgL_filled_keeling= 1/DC_Q$CO2_mgL_filled
DC_Q$CO2_mgL_filled_keeling = ifelse(is.infinite (DC_Q$CO2_mgL_filled_keeling), NA, DC_Q$CO2_mgL_filled_keeling)


ggplot(
  DC_Q %>% pivot_longer(
                     cols = c(CO2_14C_Modern, d13C_CO2),
                     names_to = "Carbon_isotope",
                     values_to = "isotope_value"
                   ),
       aes(y=isotope_value, x=CO2_mgL_filled_keeling, 
           color=Treatment))+
  
    geom_smooth(method="glm", se=F, 
                show.legend = F)+ 
    facet_wrap(~Carbon_isotope, scale="free", nrow=2)+
    stat_regline_equation(
    label.y.npc = "top", label.x.npc = "left", #0.5
    aes(label =  paste(..eq.label.., ..adj.rr.label.., sep = "~~~~")), show.legend = F, size=4)+
    
  geom_point() +
  scale_color_manual(values = treatments_colors) +
  scale_fill_manual(values = treatments_colors) +
  labs(x= "1/CO2", y="isotopic ratio")+
  theme(legend.position = "right")

Keeling plot of CO2 isotope ratio, stable and radiogenic, coloured by treatment

3.4.0.1 ANCOVA test

Does the keeling relationship for d13C-CO2 varies significantly between treatment?

👈 see code here 
#Run ANCOVA test - Difference in hydro. response between Sites
knitr::kable(DC_Q %>% anova_test(d13C_CO2 ~ CO2_mgL_filled_keeling*Treatment),
             caption="ANCOVA test results, 13C-CO2 Keeling | Treatment")
ANCOVA test results, 13C-CO2 Keeling | Treatment
Effect DFn DFd F p p<.05 ges
CO2_mgL_filled_keeling 1 36 16.913 0.000217 * 0.320
Treatment 2 36 0.154 0.858000 0.008
CO2_mgL_filled_keeling:Treatment 2 36 0.220 0.804000 0.012
👈 see code here 
#DC_Q %>% anova_test(d13C_CO2 ~ CO2_mgL_filled_keeling*Site_id)

Does the keeling relationship for 14C-CO2 varies significantly between treatment?

👈 see code here 
#Run ANCOVA test - Difference in hydro. response between Sites
knitr::kable(DC_Q %>% anova_test(CO2_14C_Modern ~ CO2_mgL_filled_keeling*Treatment),
             caption="ANCOVA test results, 14C-CO2 Keeling | Treatment")
ANCOVA test results, 14C-CO2 Keeling | Treatment
Effect DFn DFd F p p<.05 ges
CO2_mgL_filled_keeling 1 36 5.960 0.020 * 0.142
Treatment 2 36 1.821 0.176 0.092
CO2_mgL_filled_keeling:Treatment 2 36 0.308 0.737 0.017
👈 see code here 
#DC_Q %>% anova_test(CO2_14C_Modern ~ CO2_mgL_filled_keeling*Site_id)
Interpretation

  • The keeling plot suggests that CO2 concentration has a significant effect on both d13C-value and 14C-concent, which supports the hypothesis of a biological control with rapid turnover.

  • This relationship doesn’t change signficantly with treatment. Perhaps the sample size is not enough to identify a meaningful effect?

  • The d13C source of CO2 is (-25‰) (no signficant effect of slope or intercept across sites or treatment)

  • The 14C source of CO2 is between 99, 100m and 94. Which are large differences in 14C terms, but don’t appear as significant effects in the model.

References

Campeau, A., K. Bishop, N. Amvrosiadi, M. F. Billett, M. H. Garnett, H. Laudon, M. G. Öquist, and M. B. Wallin. 2019. “Current Forest Carbon Fixation Fuels Stream CO2 Emissions.” Nature Communications 10 (1). https://doi.org/10.1038/s41467-019-09922-3.
Campeau, Audrey, Kevin H. Bishop, Michael F. Billett, Mark H. Garnett, Hjalmar Laudon, Jason A. Leach, Mats B. Nilsson, Mats G. Öquist, and Marcus B. Wallin. 2017. “Aquatic Export of Young Dissolved and Gaseous Carbon from a Pristine Boreal Fen: Implications for Peat Carbon Stock Stability.” Global Change Biology 23 (12): 5523–36. https://doi.org/10.1111/gcb.13815.
Ledesma, J. L. J., T. Grabs, M. N. Futter, K. H. Bishop, H. Laudon, and S. J. Köhler. 2013. “Riparian Zone Control on Base Cation Concentration in Boreal Streams.” Biogeosciences 10 (6): 3849–68. https://doi.org/10.5194/bg-10-3849-2013.
Öquist, M. G., K. Bishop, A. Grelle, L. Klemedtsson, S. J. Köhler, H. Laudon, A. Lindroth, M. Ottosson Löfvenius, M. B. Wallin, and M. B. Nilsson. 2014. “The Full Annual Carbon Balance of Boreal Forests Is Highly Sensitive to Precipitation.” Environmental Science & Technology Letters 1 (7): 315–19. https://doi.org/10.1021/ez500169j.