Lab6_CarbonFlux
Elizabeth Olson
2026-02-15
Lab 6 Carbon Flux
In today’s lab we will be working with carboon flux data from the tower we visited last week at Oak Ridge National Lab (https://www.neonscience.org/field-sites/ornl). This lab was developed in part with resources from the NEON site https://www.neonscience.org/resources/learning-hub/tutorials/eddy-data-intro. In today’s lab you will learn how to process and interpret carbon flux data from eddie covariance towers.
Lab Submission
At the end of lab today you will upload your assignment as both a R Notebook file (.Rmd) and html file to our Canvas site. There is a template on the Canvas site to get you started with some tips on formatting you document. Your lab must be nicely organized and have a title and a header. - name your file using the format LabName_Lastname - for example for todays lab my file would be names “Lab5_Olson”
A very helpful guide on formatting R Markdown documents for assignments is available at http://www.stat.cmu.edu/~cshalizi/rmarkdown/.
Log into Posit
Go to https://docs.posit.co/cloud/get_started/. Log into your account and create a new R Studio project.
Part 1 Getting started with NEON data
First we will install and load the packages we will use. This wll take five minutes, use this tiime to read through the rest of the lab.
options(repos = c(CRAN = "https://cloud.r-project.org"))
install.packages("BiocManager")
BiocManager::install("rhdf5", ask = FALSE, update = FALSE)
install.packages("neonUtilities")
options(stringsAsFactors=F)Now lets download data from two sites the Oak Rdge Site and the Harvard Forest Site. This will take a few minutes. We will use the zipsByProduct() function from the neonUtilities package to download flux data from two sites afor one year.
knitr::opts_chunk$set(
echo = TRUE,
message = FALSE,
warning = FALSE
)
library("neonUtilities")
library("rhdf5")
zipsByProduct(dpID="DP4.00200.001", package="basic",
site=c("ORNL", "HARV"),
startdate="2024-06", enddate="2025-06",
savepath="~/Downloads",
check.size=F)## | | | 0% | |=== | 4% | |===== | 8% | |======== | 12% | |=========== | 15% | |============= | 19% | |================ | 23% | |=================== | 27% | |====================== | 31% | |======================== | 35% | |=========================== | 38% | |============================== | 42% | |================================ | 46% | |=================================== | 50% | |====================================== | 54% | |======================================== | 58% | |=========================================== | 62% | |============================================== | 65% | |================================================ | 69% | |=================================================== | 73% | |====================================================== | 77% | |========================================================= | 81% | |=========================================================== | 85% | |============================================================== | 88% | |================================================================= | 92% | |=================================================================== | 96% | |======================================================================| 100%## | | | 0% | |=== | 4% | |====== | 8% | |======== | 12% | |=========== | 16% | |============== | 20% | |================= | 24% | |==================== | 28% | |====================== | 32% | |========================= | 36% | |============================ | 40% | |=============================== | 44% | |================================== | 48% | |==================================== | 52% | |======================================= | 56% | |========================================== | 60% | |============================================= | 64% | |================================================ | 68% | |================================================== | 72% | |===================================================== | 76% | |======================================================== | 80% | |=========================================================== | 84% | |============================================================== | 88% | |================================================================ | 92% | |=================================================================== | 96% | |======================================================================| 100%There are five levels of data contained in the eddy flux bundle. For full details, refer to the NEON algorithm document.
Briefly, the data levels are:
Level 0’ (dp0p): Calibrated raw observations Level 1 (dp01): Time-aggregated observations, e.g. 30-minute mean gas concentrations Level 2 (dp02): Time-interpolated data, e.g. rate of change of a gas concentration Level 3 (dp03): Spatially interpolated data, i.e. vertical profiles Level 4 (dp04): Fluxes
We will now extract the Level 4 flux data from the HDF5 files and merge them into a single table, we’ll use the stackEddy() function from the neonUtilities package.
stackEddy() requires two inputs:
filepath: Path to a file or folder, which can be any one of: A zip file of eddy flux data downloaded from the NEON data portal A folder of eddy flux data downloaded by the zipsByProduct() function The folder of files resulting from unzipping either of 1 or 2 One or more HDF5 files of NEON eddy flux data level: dp01-4 Input the filepath you downloaded to using zipsByProduct() earlier, including the filestoStack00200 folder created by the function, and dp04:
flux <- stackEddy(filepath="~/Downloads/filesToStack00200",
level="dp04")## | | | 0% | |=== | 4% | |===== | 8% | |======== | 12% | |=========== | 15% | |============= | 19% | |================ | 23% | |=================== | 27% | |====================== | 31% | |======================== | 35% | |=========================== | 38% | |============================== | 42% | |================================ | 46% | |=================================== | 50% | |====================================== | 54% | |======================================== | 58% | |=========================================== | 62% | |============================================== | 65% | |================================================ | 69% | |=================================================== | 73% | |====================================================== | 77% | |========================================================= | 81% | |=========================================================== | 85% | |============================================================== | 88% | |================================================================= | 92% | |=================================================================== | 96% | |======================================================================| 100%## | | | 0% | |=== | 4% | |===== | 8% | |======== | 12% | |=========== | 15% | |============= | 19% | |================ | 23% | |=================== | 27% | |====================== | 31% | |======================== | 35% | |=========================== | 38% | |============================== | 42% | |================================ | 46% | |=================================== | 50% | |====================================== | 54% | |======================================== | 58% | |=========================================== | 62% | |============================================== | 65% | |================================================ | 69% | |=================================================== | 73% | |====================================================== | 77% | |========================================================= | 81% | |=========================================================== | 85% | |============================================================== | 88% | |================================================================= | 92% | |=================================================================== | 96% | |======================================================================| 100%## | | | 0% | |=================================== | 50% | |======================================================================| 100%## | | | 0% | |=================================== | 50% | |==================================================== | 75% | |======================================================================| 100%We now have extracted an object called flux. It’s a named list containing four tables: one table for each site’s data, and variables and objDesc tables. We can look at the contents of the file by calling the names() and head() function.
names(flux)## [1] "HARV" "ORNL"
## [3] "variables" "objDesc"
## [5] "issueLog" "citation_00200_RELEASE-2026"head(flux$ORNL)## timeBgn timeEnd data.fluxCo2.nsae.flux
## 1 2024-06-01 00:00:00 2024-06-01 00:29:59 NaN
## 2 2024-06-01 00:30:00 2024-06-01 00:59:59 3.136724
## 3 2024-06-01 01:00:00 2024-06-01 01:29:59 4.224824
## 4 2024-06-01 01:30:00 2024-06-01 01:59:59 7.239727
## 5 2024-06-01 02:00:00 2024-06-01 02:29:59 6.542591
## 6 2024-06-01 02:30:00 2024-06-01 02:59:59 4.640835
## data.fluxCo2.stor.flux data.fluxCo2.turb.flux data.fluxH2o.nsae.flux
## 1 NaN -1.394239 NaN
## 2 1.5821981 1.554526 12.387344
## 3 0.9571699 3.267654 11.354807
## 4 4.0928084 3.146919 28.375953
## 5 2.0369123 4.505679 27.026174
## 6 1.3765265 3.264309 6.830413
## data.fluxH2o.stor.flux data.fluxH2o.turb.flux data.fluxMome.turb.veloFric
## 1 NaN 36.820222 0.3509450
## 2 -22.645517 35.032861 0.3290410
## 3 -3.128678 14.483484 0.3268819
## 4 18.574514 9.801439 0.3452148
## 5 -6.302745 33.328920 0.4300269
## 6 -6.579157 13.409570 0.3271710
## data.fluxTemp.nsae.flux data.fluxTemp.stor.flux data.fluxTemp.turb.flux
## 1 -47.34805 -9.306308 -38.04174
## 2 -60.65596 -3.493420 -57.16254
## 3 -59.07815 -9.831581 -49.24657
## 4 -62.55373 -10.891305 -51.66243
## 5 -93.91015 -3.250511 -90.65963
## 6 -49.41696 2.575467 -51.99243
## data.foot.stat.angZaxsErth data.foot.stat.distReso
## 1 52.47689 20
## 2 58.40503 20
## 3 62.01288 20
## 4 68.81185 20
## 5 73.06157 20
## 6 76.44452 20
## data.foot.stat.veloYaxsHorSd data.foot.stat.veloZaxsHorSd
## 1 0.6314968 0.4512513
## 2 0.6276120 0.4547496
## 3 0.6019615 0.4234288
## 4 0.6609377 0.4970567
## 5 0.7172832 0.5274381
## 6 0.5831353 0.4506983
## data.foot.stat.veloFric data.foot.stat.distZaxsMeasDisp
## 1 0.3507455 21.6885
## 2 0.3290410 21.6885
## 3 0.3266543 21.6885
## 4 0.3452148 21.6885
## 5 0.4300269 21.6885
## 6 0.3268820 21.6885
## data.foot.stat.distZaxsRgh data.foot.stat.distObkv data.foot.stat.paraStbl
## 1 1.581997 106.67445 0.2033148
## 2 2.168850 56.88289 0.3812834
## 3 2.168850 63.00049 0.3442592
## 4 2.168850 70.27310 0.3086316
## 5 2.168850 78.75587 0.2753890
## 6 2.168850 59.66478 0.3635059
## data.foot.stat.distZaxsAbl data.foot.stat.distXaxs90
## 1 1000 420
## 2 1000 360
## 3 1000 360
## 4 1000 340
## 5 1000 380
## 6 1000 360
## data.foot.stat.distXaxsMax data.foot.stat.distYaxs90 qfqm.fluxCo2.nsae.qfFinl
## 1 180 100 1
## 2 160 100 0
## 3 160 100 0
## 4 160 100 1
## 5 160 100 1
## 6 160 100 0
## qfqm.fluxCo2.stor.qfFinl qfqm.fluxCo2.turb.qfFinl qfqm.fluxH2o.nsae.qfFinl
## 1 1 0 1
## 2 0 0 0
## 3 0 0 0
## 4 1 0 1
## 5 1 0 1
## 6 0 0 0
## qfqm.fluxH2o.stor.qfFinl qfqm.fluxH2o.turb.qfFinl qfqm.fluxMome.turb.qfFinl
## 1 1 0 0
## 2 0 0 0
## 3 0 0 0
## 4 1 0 0
## 5 1 0 0
## 6 0 0 0
## qfqm.fluxTemp.nsae.qfFinl qfqm.fluxTemp.stor.qfFinl qfqm.fluxTemp.turb.qfFinl
## 1 0 0 0
## 2 0 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 0 0
## qfqm.foot.turb.qfFinl
## 1 0
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0We can see the units and variables by calling the term ‘variables’
flux$variables## category system variable stat units
## 1 data fluxCo2 nsae timeBgn NA
## 2 data fluxCo2 nsae timeEnd NA
## 3 data fluxCo2 nsae flux umolCo2 m-2 s-1
## 4 data fluxCo2 stor timeBgn NA
## 5 data fluxCo2 stor timeEnd NA
## 6 data fluxCo2 stor flux umolCo2 m-2 s-1
## 7 data fluxCo2 turb timeBgn NA
## 8 data fluxCo2 turb timeEnd NA
## 9 data fluxCo2 turb flux umolCo2 m-2 s-1
## 10 data fluxH2o nsae timeBgn NA
## 11 data fluxH2o nsae timeEnd NA
## 12 data fluxH2o nsae flux W m-2
## 13 data fluxH2o stor timeBgn NA
## 14 data fluxH2o stor timeEnd NA
## 15 data fluxH2o stor flux W m-2
## 16 data fluxH2o turb timeBgn NA
## 17 data fluxH2o turb timeEnd NA
## 18 data fluxH2o turb flux W m-2
## 19 data fluxMome turb timeBgn NA
## 20 data fluxMome turb timeEnd NA
## 21 data fluxMome turb veloFric m s-1
## 22 data fluxTemp nsae timeBgn NA
## 23 data fluxTemp nsae timeEnd NA
## 24 data fluxTemp nsae flux W m-2
## 25 data fluxTemp stor timeBgn NA
## 26 data fluxTemp stor timeEnd NA
## 27 data fluxTemp stor flux W m-2
## 28 data fluxTemp turb timeBgn NA
## 29 data fluxTemp turb timeEnd NA
## 30 data fluxTemp turb flux W m-2
## 31 data foot stat timeBgn NA
## 32 data foot stat timeEnd NA
## 33 data foot stat angZaxsErth deg
## 34 data foot stat distReso m
## 35 data foot stat veloYaxsHorSd m s-1
## 36 data foot stat veloZaxsHorSd m s-1
## 37 data foot stat veloFric m s-1
## 38 data foot stat distZaxsMeasDisp m
## 39 data foot stat distZaxsRgh m
## 40 data foot stat distObkv m
## 41 data foot stat paraStbl -
## 42 data foot stat distZaxsAbl m
## 43 data foot stat distXaxs90 m
## 44 data foot stat distXaxsMax m
## 45 data foot stat distYaxs90 m
## 46 qfqm fluxCo2 nsae timeBgn NA
## 47 qfqm fluxCo2 nsae timeEnd NA
## 48 qfqm fluxCo2 nsae qfFinl NA
## 49 qfqm fluxCo2 stor qfFinl NA
## 50 qfqm fluxCo2 stor timeBgn NA
## 51 qfqm fluxCo2 stor timeEnd NA
## 52 qfqm fluxCo2 turb timeBgn NA
## 53 qfqm fluxCo2 turb timeEnd NA
## 54 qfqm fluxCo2 turb qfFinl NA
## 55 qfqm fluxH2o nsae timeBgn NA
## 56 qfqm fluxH2o nsae timeEnd NA
## 57 qfqm fluxH2o nsae qfFinl NA
## 58 qfqm fluxH2o stor qfFinl NA
## 59 qfqm fluxH2o stor timeBgn NA
## 60 qfqm fluxH2o stor timeEnd NA
## 61 qfqm fluxH2o turb timeBgn NA
## 62 qfqm fluxH2o turb timeEnd NA
## 63 qfqm fluxH2o turb qfFinl NA
## 64 qfqm fluxMome turb timeBgn NA
## 65 qfqm fluxMome turb timeEnd NA
## 66 qfqm fluxMome turb qfFinl NA
## 67 qfqm fluxTemp nsae timeBgn NA
## 68 qfqm fluxTemp nsae timeEnd NA
## 69 qfqm fluxTemp nsae qfFinl NA
## 70 qfqm fluxTemp stor qfFinl NA
## 71 qfqm fluxTemp stor timeBgn NA
## 72 qfqm fluxTemp stor timeEnd NA
## 73 qfqm fluxTemp turb timeBgn NA
## 74 qfqm fluxTemp turb timeEnd NA
## 75 qfqm fluxTemp turb qfFinl NA
## 76 qfqm foot turb timeBgn NA
## 77 qfqm foot turb timeEnd NA
## 78 qfqm foot turb qfFinl NAQuestion 1
What unit is the carbon dioxide flux in within the dataset?
## Question 2
What time of year is carbon flux variability lowest? Why do you think that is?
Now lets look at daily cycling.
# First we want to make sure time is POSIXct
flux$ORNL$timeBgn <- as.POSIXct(flux$ORNL$timeBgn, tz = "GMT")
# Then we extract components
flux$ORNL$year <- format(flux$ORNL$timeBgn, "%Y")
flux$ORNL$yday <- format(flux$ORNL$timeBgn, "%Y-%m-%d")
flux$ORNL$hour <- as.numeric(format(flux$ORNL$timeBgn, "%H"))
#average by day
daily_mean <- aggregate(
data.fluxCo2.nsae.flux ~ yday,
data = flux$ORNL,
FUN = mean,
na.rm = TRUE
)
#average by hour
hourly_year_mean <- aggregate(
data.fluxCo2.nsae.flux ~ year + hour,
data = flux$ORNL,
FUN = mean,
na.rm = TRUE
)
#plot yearly cycle
# Get unique years
yrs <- unique(hourly_year_mean$year)
# Set up plot
plot(NULL,
xlim = c(0, 23),
ylim = range(hourly_year_mean$data.fluxCo2.nsae.flux, na.rm = TRUE),
xlab = "Hour of Day",
ylab = "Mean CO2 Flux")
# Add each year as a line
cols <- rainbow(length(yrs))
for (i in seq_along(yrs)) {
yr_data <- subset(hourly_year_mean, year == yrs[i])
lines(yr_data$hour,
yr_data$data.fluxCo2.nsae.flux,
col = cols[i],
lwd = 2)
}
legend("topright", legend = yrs, col = cols, lwd = 2)
Question 3
How does carbon flux vary over the day? Form a hypothesis as to what is driving this diurnal cycling? Recall that at NEON sites all data are reported using UTC time. Or Universal Time Central not University of Tennessee Chattanooga time :). Keep in mind that to convert US-EST to UTC you add five hours. Keep the time difference in mind when interpreting the plot above.
You may have mentioned photosynthetic activity as a major driver of the cycling we see above. At NEON sites PAR or Photosynthetically Active Radiation is measured as well and we can compare our flux data to these data.
#download PAR data
pr <- loadByProduct("DP1.00024.001", site="ORNL",
timeIndex=30, package="basic",
startdate="2024-06", enddate="2025-06",
check.size=F)## | | | 0% | |===== | 8% | |=========== | 15% | |================ | 23% | |====================== | 31% | |=========================== | 38% | |================================ | 46% | |====================================== | 54% | |=========================================== | 62% | |================================================ | 69% | |====================================================== | 77% | |=========================================================== | 85% | |================================================================= | 92% | |======================================================================| 100%## | | | 0% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |======= | 9% | |======= | 10% | |======== | 11% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |================ | 22% | |================ | 23% | |================= | 24% | |================== | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%## | | | 0% | |======================= | 33% | |=============================================== | 67% | |======================================================================| 100%Now lets link both data sets to the same time step.
#connect PAR data from the tower top to the flux data
pr.top <- pr$PARPAR_30min[which(pr$PARPAR_30min$verticalPosition==
max(pr$PARPAR_30min$verticalPosition)),]
#create new variable
pr.top$timeBgn <- pr.top$startDateTime
#merge data
fx.pr <- merge(pr.top, flux$ORNL, by="timeBgn")Now we can plot carbon flux versus PAR, to test our hypothesis.
plot(fx.pr$data.fluxCo2.nsae.flux~fx.pr$PARMean,
pch=".", ylim=c(-50,50),
xlab="PAR", ylab="CO2 flux")
Question 4
Do PAR and carbon flux correlate? How might you test this in a more quantitative way?
Lets look at the least regression between the two data sets.
# Define variables
ts2 <- fx.pr$data.fluxCo2.nsae.flux
ts1 <- fx.pr$PARMean
# Remove NA pairs
df_lm <- data.frame(ts1, ts2)
df_lm <- df_lm[complete.cases(df_lm), ]
# Run linear model
mod <- lm(ts2 ~ ts1, data = df_lm)
# Extract statistics
r <- cor(df_lm$ts1, df_lm$ts2)
r2 <- summary(mod)$r.squared
p <- summary(mod)$coefficients[2, 4]
# Print results
cat("r =", round(r, 3), "\n")## r = -0.416cat("R² =", round(r2, 3), "\n")## R² = 0.173cat("p-value =", signif(p, 3), "\n")## p-value = 0Question 5
Is carbon flux (CO2) significantly correlated with PAR ? What other factors might influence this correlation?
Part 2 Vertical Profiles
NEON sites also have spatially interpolated profiles from the towers as well as the rates of change of H2O, and isotope values. Lets look at these variables and their relationship with CO2 flux. This will take a few minutes.
iso <- stackEddy(filepath="~/Downloads/filesToStack00200/",
level="dp01", var=c("rtioMoleDryCo2","rtioMoleDryH2o",
"dlta13CCo2","dlta18OH2o"), avg=30)## | | | 0% | |=== | 4% | |===== | 8% | |======== | 12% | |=========== | 15% | |============= | 19% | |================ | 23% | |=================== | 27% | |====================== | 31% | |======================== | 35% | |=========================== | 38% | |============================== | 42% | |================================ | 46% | |=================================== | 50% | |====================================== | 54% | |======================================== | 58% | |=========================================== | 62% | |============================================== | 65% | |================================================ | 69% | |=================================================== | 73% | |====================================================== | 77% | |========================================================= | 81% | |=========================================================== | 85% | |============================================================== | 88% | |================================================================= | 92% | |=================================================================== | 96% | |======================================================================| 100%## | | | 0% | |=== | 4% | |===== | 8% | |======== | 12% | |=========== | 15% | |============= | 19% | |================ | 23% | |=================== | 27% | |====================== | 31% | |======================== | 35% | |=========================== | 38% | |============================== | 42% | |================================ | 46% | |=================================== | 50% | |====================================== | 54% | |======================================== | 58% | |=========================================== | 62% | |============================================== | 65% | |================================================ | 69% | |=================================================== | 73% | |====================================================== | 77% | |========================================================= | 81% | |=========================================================== | 85% | |============================================================== | 88% | |================================================================= | 92% | |=================================================================== | 96% | |======================================================================| 100%## | | | 0% | |=================================== | 50% | |======================================================================| 100%## | | | 0% | |=================================== | 50% | |==================================================== | 75% | |======================================================================| 100%Next we will use the time stamps in a different way, using aggregate() to select all of the records for a single year. And discard the verticalPosition values that are string values - those are the calibration gases.
iso.d <- iso$ORNL # use Oak Ridge dataset
# Clean height column
iso.d$verticalPosition <- as.numeric(iso.d$verticalPosition)
iso.d <- iso.d[!is.na(iso.d$verticalPosition), ]
#combine data using aggregate function
profile.avg <- aggregate(
data.co2Stor.rtioMoleDryCo2.mean ~ verticalPosition,
data = iso.d,
FUN = function(x) c(mean = mean(x, na.rm = TRUE),
sd = sd(x, na.rm = TRUE))
)
#make data matrix for plotting
profile.avg <- do.call(data.frame, profile.avg)
names(profile.avg) <- c("verticalPosition", "mean_co2", "sd_co2")
profile.avg <- profile.avg[order(profile.avg$verticalPosition), ]
#plot data
plot(profile.avg$mean_co2,
profile.avg$verticalPosition,
type = "b",
pch = 19,
xlab = "Mean CO2",
ylab = "Tower level")
arrows(profile.avg$mean_co2 - profile.avg$sd_co2,
profile.avg$verticalPosition,
profile.avg$mean_co2 + profile.avg$sd_co2,
profile.avg$verticalPosition,
angle = 90,
code = 3,
length = 0.05)
## Question 6
How does the mean CO2 concentration change with height about the ground? What causes this?
Now lets calculate the vertical profile for the Harvard forest site and compare them.
iso.d <- iso$HARV # use Harvard dataset
# Clean height column
iso.d$verticalPosition <- as.numeric(iso.d$verticalPosition)
iso.d <- iso.d[!is.na(iso.d$verticalPosition), ]
#combine data using aggregate function
profile.avg2 <- aggregate(
data.co2Stor.rtioMoleDryCo2.mean ~ verticalPosition,
data = iso.d,
FUN = function(x) c(mean = mean(x, na.rm = TRUE),
sd = sd(x, na.rm = TRUE))
)
#make data matrix for plotting
profile.avg2 <- do.call(data.frame, profile.avg2)
names(profile.avg2) <- c("verticalPosition", "mean_co2", "sd_co2")
profile.avg2 <- profile.avg2[order(profile.avg2$verticalPosition), ]
# Make plots. First determine axis limits
xrange <- range(
c(profile.avg$mean_co2 + profile.avg$sd_co2,
profile.avg$mean_co2 - profile.avg$sd_co2,
profile.avg2$mean_co2 + profile.avg2$sd_co2,
profile.avg2$mean_co2 - profile.avg2$sd_co2),
na.rm = TRUE
)
yrange <- range(
c(profile.avg$verticalPosition,
profile.avg2$verticalPosition)
)
plot(NULL,
xlim = xrange,
ylim = yrange,
xlab = "Mean CO2",
ylab = "Tower level",
main = "Average Vertical CO2 Profile")
# ORNL (red)
lines(profile.avg$mean_co2,
profile.avg$verticalPosition,
type = "b",
pch = 19,
col = "red",
lwd = 2)
arrows(profile.avg$mean_co2 - profile.avg$sd_co2,
profile.avg$verticalPosition,
profile.avg$mean_co2 + profile.avg$sd_co2,
profile.avg$verticalPosition,
angle = 90, code = 3, length = 0.05,
col = "red")
# HARV (blue)
lines(profile.avg2$mean_co2,
profile.avg2$verticalPosition,
type = "b",
pch = 19,
col = "blue",
lwd = 2)
arrows(profile.avg2$mean_co2 - profile.avg2$sd_co2,
profile.avg2$verticalPosition,
profile.avg2$mean_co2 + profile.avg2$sd_co2,
profile.avg2$verticalPosition,
angle = 90, code = 3, length = 0.05,
col = "blue")
legend("topright",
legend = c("ORNL", "HARV"),
col = c("red", "blue"),
pch = 19,
lwd = 2,
bty = "n")
Now lets plot the data.
Question 7
How does carbon flux vary between the two sites? What is a potential reason for this difference?
Now lets compare the ORNL and HARV sites carbon isotope values
iso.d <- iso$ORNL # use Oak Ridge dataset
# Clean height column
iso.d$verticalPosition <- as.numeric(iso.d$verticalPosition)
iso.d <- iso.d[!is.na(iso.d$verticalPosition), ]
#combine data using aggregate function
profile.avg <- aggregate(
data.isoCo2.dlta13CCo2.mean ~ verticalPosition,
data = iso.d,
FUN = function(x) c(mean = mean(x, na.rm = TRUE),
sd = sd(x, na.rm = TRUE))
)
#make data matrix for plotting
profile.avg <- do.call(data.frame, profile.avg)
names(profile.avg) <- c("verticalPosition", "mean_d13c", "sd_d13c")
profile.avg <- profile.avg[order(profile.avg$verticalPosition), ]
iso.d <- iso$HARV # use Harvard dataset
# Clean height column
iso.d$verticalPosition <- as.numeric(iso.d$verticalPosition)
iso.d <- iso.d[!is.na(iso.d$verticalPosition), ]
#combine data using aggregate function
profile.avg2 <- aggregate(
data.isoCo2.dlta13CCo2.mean ~ verticalPosition,
data = iso.d,
FUN = function(x) c(mean = mean(x, na.rm = TRUE),
sd = sd(x, na.rm = TRUE))
)
#make data matrix for plotting
profile.avg2 <- do.call(data.frame, profile.avg2)
names(profile.avg2) <- c("verticalPosition", "mean_d13c", "sd_d13c")
profile.avg2 <- profile.avg2[order(profile.avg2$verticalPosition), ]
# Make plots. First determine axis limits
xrange <- range(
c(profile.avg$mean_d13c + profile.avg$sd_d13c,
profile.avg$mean_d13c - profile.avg$sd_d13c,
profile.avg2$mean_d13c + profile.avg2$sd_d13c,
profile.avg2$mean_d13c - profile.avg2$sd_d13c),
na.rm = TRUE
)
yrange <- range(
c(profile.avg$verticalPosition,
profile.avg2$verticalPosition)
)
plot(NULL,
xlim = xrange,
ylim = yrange,
xlab = "Mean d13c",
ylab = "Tower level",
main = "Average Vertical d13c Profile")
# ORNL (red)
lines(profile.avg$mean_d13c,
profile.avg$verticalPosition,
type = "b",
pch = 19,
col = "maroon",
lwd = 2)
arrows(profile.avg$mean_d13c - profile.avg$sd_d13c,
profile.avg$verticalPosition,
profile.avg$mean_d13c + profile.avg$sd_d13c,
profile.avg$verticalPosition,
angle = 90, code = 3, length = 0.05,
col = "maroon")
# HARV (blue)
lines(profile.avg2$mean_d13c,
profile.avg2$verticalPosition,
type = "b",
pch = 19,
col = "lightblue",
lwd = 2)
arrows(profile.avg2$mean_d13c - profile.avg2$sd_d13c,
profile.avg2$verticalPosition,
profile.avg2$mean_d13c + profile.avg2$sd_d13c,
profile.avg2$verticalPosition,
angle = 90, code = 3, length = 0.05,
col = "lightblue")
legend("topright",
legend = c("ORNL", "HARV"),
col = c("maroon", "lightblue"),
pch = 19,
lwd = 2,
bty = "n")
Question 8
How do the average carbon-isotope delta values differ between the two sites? What about these two sites might cause this?
Question 9
Graduate level (undergrad optional) Pick another site from the NEON network (https://www.neonscience.org/field-sites) to compare data from. How do the carbon isotope values differ between these sites, why are these differences larger or smaller than the comparison above in question 8?
# include your code here for your answer to queston 9