Fluxes

Climatic data for Paracou and Tapajos used to calibrate the new leaf phenology module in TROLL version 4 and evaluate simulation results.

Carbon and water fluxes - FLUXNET

FLUXNET 2015 (Pastorello, Gilberto et al. 2020) at Paracou and Tapajos was retrieved freely from FLUXNET (the manifest is available in data). They are mainly used for TROLL model forcing but also include fluxes for TROLL simulations validation including growth primary productivity (GPP) and evapotranspiration (ET). We also computed Photosynthetically Active Radiation (PAR) as it is a driver of GPP, and daily maximum Vapour Pressure Deficit (VPD) as it is a driver of ET. We removed GPP and ET data between 2006 and 2009 in Tapajos due to a known incident with the sensors.

The FLUXNET2015 Dataset includes data collected at sites from multiple regional flux networks. The preparation of this FLUXNET Dataset has been possible thanks only to the efforts of many scientists and technicians around the world and the coordination among teams from regional networks.

PI: Damien Bonal (Paracou), Mike Goulden (Tapajos), Scott Saleska (Tapajos)

Code 
lambda <- function(t_air) 2.501 - (2.361*10^-3)*t_air
fluxnet <- list(
  "Paracou" = "data/climate/FLX_GF-Guy_FLUXNET2015_SUBSET_HH_2004-2014_2-4.csv",
  "Tapajos" = "data/climate/FLX_BR-Sa1_FLUXNET2015_SUBSET_HR_2002-2011_1-4.csv"
) %>% 
  lapply(vroom::vroom, na = "-9999") %>% 
  bind_rows(.id = "site") %>% 
  select(site, TIMESTAMP_START, TIMESTAMP_END, P_F, SW_IN_F, 
         TA_F, VPD_F, WS_F, GPP_NT_VUT_REF, LE_F_MDS) %>% 
  mutate(time = as_datetime(as.character(as.POSIXlt(as.character(TIMESTAMP_START), 
                                                    format = "%Y%m%d%H%M")))) %>% 
  rename(rainfall = P_F, snet = SW_IN_F, temperature = TA_F, vpd = VPD_F, ws = WS_F,
         gpp = GPP_NT_VUT_REF, le = LE_F_MDS) %>% 
  mutate(et = 10^-6*le*60*30 / lambda(temperature)) %>% # W.m-2.half-hour to MJ.s-1 to kg.m-2.half-hour
  mutate(et = ifelse(hour(time) %in% 7:18, et, NA)) %>% 
  mutate(vpd = vpd/10,
         gpp = gpp/100, par = snet*2.27*60*30/10^3) %>% 
  select(site, time, vpd, par, gpp, et, rainfall, ws,temperature)
fluxnet_tapajos <- filter(fluxnet, site == "Tapajos")
era <- vroom::vroom("data/climate/era5_tapajos.csv") %>% 
  filter(year(date) > 2001) %>% 
  mutate(time = as_datetime(date), rainfall = precipitation) %>% 
  select(time, rainfall)
fluxnet_tapajos <- fluxnet_tapajos %>% 
  select(-rainfall) %>% 
  left_join(era)
time_freq <- 0.5 # TROLL v4
start <- as_datetime(as_date(min(fluxnet_tapajos$time))) # to start from 00:30
stop <- max(fluxnet_tapajos$time) + 30*60 # to stop at 23:00
time <- seq(start, stop, by = 60 * 60 * time_freq)
data <- left_join(tibble(time = time),
    fluxnet_tapajos,
    by = "time"
  ) %>%
  group_by(date = as_date(time)) %>%
  mutate(across(c(temperature, vpd, ws, gpp),
                ~ zoo::na.spline(., time, na.rm = FALSE))) %>% 
  mutate(across(c(par, et),
                ~ zoo::na.approx(., time, na.rm = FALSE))) %>% 
  ungroup() %>% 
  select(-date) %>% 
  mutate(across(c(par), ~ ifelse(is.na(.), 0, .))) %>% 
  mutate(across(c(par), ~ ifelse(. < 0, 0, .))) %>% 
  mutate(site = "Tapajos")
bind_rows(filter(fluxnet, site == "Paracou"), data) %>% 
    mutate(vpd_max = vpd) %>% 
  select(-vpd) %>% 
  gather(variable, value, -site, -time) %>% 
  group_by(site, variable, date = date(time)) %>%
  select(-time) %>% 
  mutate(value = ifelse(value %in% c("et", "par", "rainfall"),
                           sum(value, na.rm = T), value)) %>% 
  mutate(value = ifelse(variable %in% c("vpd_max"),
                           max(value, na.rm = T), value)) %>% 
  summarise_all(mean, na.rm = T) %>% 
  mutate(value = ifelse(variable %in% c("et", "gpp") & 
                          site == "Tapajos" & year(date) %in% 2006:2009,
                           NA, value)) %>% 
  write_tsv("outputs/fluxnet_fluxes.tsv")
Code 
read_tsv("outputs/fluxnet_fluxes.tsv") %>% 
  ggplot(aes(date, value, col = site)) +
  geom_line(alpha = 0.5) +
  geom_line(data = read_tsv("outputs/fluxnet_fluxes.tsv") %>% 
  group_by(variable, site, date = floor_date(date, "month")) %>% 
  summarise_all(mean)) +
  theme_bw() +
  facet_wrap(~ variable, scales = "free_y") +
  theme(axis.title = element_blank())

Yearly means of fluxes at Paracou and Tapajos across ten years.
Code 
floor_biweek <- function(full) 
  as_date(paste0(year(full), "-", month(full), "-", ifelse(day(full) > 15, 15, 1)))
read_tsv("outputs/fluxnet_fluxes.tsv") %>% 
  group_by(site, variable, date = floor_date(date, "15 days")) %>% 
  summarise_all(mean, na.rm = T) %>%
  mutate(week = week(floor_biweek(date))) %>%  
  group_by(site, variable, week) %>% 
  summarise(l = quantile(value, 0.025, na.rm = TRUE), 
            m = mean(value, na.rm = TRUE), 
            h = quantile(value, 0.975, na.rm = TRUE)) %>% 
  ggplot(aes(week, m, col = site, fill = site)) +
  geom_ribbon(aes(ymin = l, ymax = h), col = NA, alpha = 0.2) +
  geom_smooth(se = FALSE) +
  theme_bw() +
  theme(axis.title = element_blank()) +
  ggtitle("15-days means variations") +
  scale_fill_discrete("") +
  scale_color_discrete("") +
  facet_wrap(~ variable, scales = "free_y") +
    scale_x_continuous(breaks = week(as_date(paste("2000-", 1:12, "-01"))),
                     labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"))

Mean annual cycle of fortnightly means of fluxes at Paracou and Tapajos across ten years.
Code 
list(
  "Paracou" = "data/climate/FLX_GF-Guy_FLUXNET2015_SUBSET_HH_2004-2014_2-4.csv",
  "Tapajos" = "data/climate/FLX_BR-Sa1_FLUXNET2015_SUBSET_HR_2002-2011_1-4.csv"
) %>% 
  lapply(vroom::vroom, na = "-9999") %>% 
  bind_rows(.id = "site") %>% 
  select(site, TIMESTAMP_START, TIMESTAMP_END, P_F, SW_IN_F, 
         TA_F, VPD_F, WS_F, GPP_NT_VUT_REF, LE_F_MDS) %>% 
  mutate(time = as_datetime(as.character(as.POSIXlt(as.character(TIMESTAMP_START), 
                                                    format = "%Y%m%d%H%M")))) %>% 
  rename(gpp = GPP_NT_VUT_REF) %>% 
  mutate(gpp = gpp/100*10^2*365/10^3) %>% 
  select(site, time, gpp) %>% 
  mutate(d_start = ifelse(site == "Paracou", yday("2020/08/01"), yday("2020/06/15"))) %>% 
  mutate(d_end = ifelse(site == "Paracou", yday("2020/12/01"), yday("2020/11/01"))) %>% 
  mutate(season = ifelse(yday(time) %in% d_start:d_end, "dry", "wet")) %>% 
  group_by(site, season, time = hms::as_hms(time)) %>% 
  summarise(l = quantile(gpp, 0.025, na.rm = TRUE), 
            m = mean(gpp, na.rm = TRUE), 
            h = quantile(gpp, 0.975, na.rm = TRUE)) %>% 
  mutate(time = as.numeric(time)/(60*60)) %>% 
  ggplot(aes(time, m, col = season, fill = season)) +
  geom_ribbon(aes(ymin = l, ymax = h), col = NA, alpha = 0.2) +
  geom_line() +
  geom_point() +
  theme_bw() +
  facet_wrap(~ site) +
  theme(legend.position = "bottom") +
  ylab(expression("Gross Primary Productivity ["~kgC~m^{-2}~yr^{-1}~"]")) +
  xlab("")

GPP half-hourly.

Carbon flux - TROPOMI SIF

Tapajos and Paracou Solar Induced Fluorescence (SIF) from TROPOMI satellites was retrieved from https://figshare.com/ndownloader/articles/19336346/versions/2 and processed with the R script available in data. The SIF values were converted in Growth Primary Productivity (GPP) using \(GPP = 15.343 \times SIF\) following Chen et al. (2022).

Code 
read_tsv("data/fluxes/rtsif.tsv") %>% 
  mutate(value = rtsif*15.343*365/10^3) %>% 
  ggplot(aes(date, value)) +
  geom_line(alpha = 0.5, aes(col = site)) +
  theme_bw() +
  theme(legend.position = "bottom") +
  xlab("") + ylab("Growth Primary Productivity (kgC/m2/yr)")

8-days variations of GPP derived from TROPOMI SIF at Paracou and Tapajos across twenty years.
Code 
read_tsv("data/fluxes/rtsif.tsv") %>% 
  mutate(value = rtsif*15.343*365/10^3) %>% 
  mutate(week = week(date)) %>% 
  group_by(site, week) %>% 
  summarise(l = quantile(value*5, 0.025, na.rm = TRUE), 
            m = quantile(value*5, 0.5, na.rm = TRUE), 
            h = quantile(value*5, 0.975, na.rm = TRUE)) %>% 
  ggplot(aes(week, m, col = site, fill = site)) +
  geom_ribbon(aes(ymin = l, ymax = h), col = NA, alpha = 0.2) +
  geom_smooth(se = FALSE) +
  geom_point() +geom_line() +
  theme_bw() +
  scale_fill_discrete("") +
  scale_color_discrete("") +
    scale_x_continuous(breaks = week(as_date(paste("2000-", 1:12, "-01"))),
                     labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D")) +
  xlab("") + ylab("Growth Primary Productivity (kgC/m2/yr)")

Mean annual cycle of 8-days GPP derived from TROPOMI SIF at Paracou and Tapajos across ten years.

Soil water content

Guyaflux are Paracou Eddy Flux tower raw data (BONAL et al. 2008) from 2004 to 2022. Vapour Pressure Deficit (VPD, kPa) is computed from relative humidity (%). Data are available on demand and can be partly retrieved from PLUMBER2 (Ukkola, Abramowitz, and De Kauwe 2022) or FLUXNET (Pastorello, Gilberto et al. 2020). Guyaflux was mainly used for TROLL simulations validation including soil water content (SWC) at different depths. TROLL 4.0 water fluxes were assessed using the relative variation of soil water content (RSWC, %) from the top horizon from Paracou eddy flux tower (Bonal et al., 2008) and from reanalysed top horizon eddy flux tower against climatic water deficit at Tapajos (Restrepo-Coupe et al., 2024), defined as the daily mean value of soil water content (m3.m-3) divided by the yearly 95th quantile of daily mean.

PI: Damien Bonal & Natalia Restrepo-Coupe.

Code 
swc_par <- read_csv2("data/climate/GX-METEO^MEDDY-2004-2022.csv",
      locale = locale(decimal_mark = ",")) %>% 
  group_by(Year) %>% 
  mutate(date = as_date(paste0(min(Year), "/", min(Month),
                               "/", min(`Julian Day`))) + (`Julian Day` - 1)) %>% 
    separate(`Hour/min`, c("hour", "minute")) %>% 
  ungroup() %>% 
  mutate(minute = ifelse(is.na(minute), 0, 60*as.numeric(minute)/10)) %>% 
  mutate(time = as_datetime(paste(as.character(date), 
                                  paste0(hour, ":", minute, ":00")))) %>% 
  mutate(time = time - 30*60) %>% 
  filter(year(time)  %in% 2004:2014) %>% 
  mutate(swc10 = `SWC (-10 cm)`,
         swc20 = `SWC (-20 cm)`,
         swc80 = `SWC (-80 cm)`,
         swc160 = `SWC (-160 cm)`,
         swc260 = `SWC (-260 cm)`) %>%
  mutate(site = "Paracou") %>% 
  select(site, time, swc10, swc20, swc80, swc160, swc260) %>% 
  group_by(site, date = date(time)) %>% 
  summarise_all(mean, na.rm = TRUE) %>% 
  select(-time) %>% 
  gather(variable, value, -site, -date) %>% 
  na.omit()
swc_tap <- read_csv("data/fluxes/km67.eddyflux.seasonal.txt") %>% 
  mutate(date = as_date(paste0(YYYY, "-", MM, "-", DD))) %>%
  select(date, SMshallow, SMdeep) %>% 
  rename(swc10 = SMshallow, swc160 = SMdeep) %>% 
  mutate(site = "Tapajos") %>% 
  gather(variable, value, -site, -date) %>% 
  na.omit()
bind_rows(swc_par, swc_tap) %>% 
  write_tsv("outputs/guyaflux_swc.tsv")
Code 
read_tsv("outputs/guyaflux_swc.tsv") %>% 
  mutate(depth = as.numeric(gsub("swc", "", variable))) %>% 
  group_by(site, depth, day = yday(date)) %>% 
  summarise(l = quantile(value, 0.025, na.rm = TRUE), 
            m = mean(value, na.rm = TRUE), 
            h = quantile(value, 0.975, na.rm = TRUE)) %>% 
  ggplot(aes(day, m, col = depth, fill = depth, group = depth)) +
  geom_ribbon(aes(ymin = l, ymax = h), col = NA, alpha = 0.2) +
  geom_line() +
  theme_bw() +
  scale_color_viridis_c(direction = -1) +
  scale_fill_viridis_c(direction = -1) +
  theme(axis.title = element_blank()) +
  facet_wrap(~ site) +
  ggtitle("Daily means variations of SWC") +
  scale_x_continuous(breaks = c(1, 32, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336),
                   labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"))

Mean annual cycle of daily means of soil water content.

Litterfall

Paracou litterfall were provided by Damien Bonal (D. Bonal pers. com.). Tapajos litterfall are available online through the Oak Ridge National Laboratory (ORNL) Distributed Active Archive Center (DAAC): https://daac.ornl.gov/LBA/guides/CD10_Litter_Tapajos.html . Litterfall data were mainly used for the calibration of the new leaf phenology module in TROLL version 4. Leaf litterfall data were converted in fluxes over a defined time lag (depending on site and sampling date). For calibration two indices were derived for each years:

  • Day of the litterfall pike as the day of the maximum annual value

  • The litterfall rate is the ratio between the pike value (average of the annual maximum value and the previous and following values) divided by the base flow (average values between January and April)

Code 
prop_paracou <- vroom::vroom("data/fluxes/prop_litterfall.csv", n_max = 4) %>% 
  filter(date == "feuille") %>% 
  select(-date) %>% 
  gather(date, prop_leaf)  %>% 
  mutate(date = as_date(date, format = "%d/%m/%Y")) %>% 
  na.omit() %>% 
  group_by(month = month(date)) %>% 
  summarise(prop_leaf = quantile(prop_leaf, 0.5, na.rm = TRUE)/100)
litter_paracou <- vroom::vroom("data/fluxes/litterfall_2003_2022.csv") %>% 
  rename("collector" = "Numero collecteur", "subplot" = "Carre guyaflux") %>% 
  gather(date, litter, -collector, -subplot) %>% 
  mutate(date = as_date(date, format = "%d/%m/%Y")) %>% 
  na.omit() %>% 
  group_by(date) %>% 
  summarise(litter = mean(litter)) %>% 
  mutate(dt = as.numeric(date - lag(date))) %>% 
  mutate(month = month(date)) %>% 
  left_join(prop_paracou) %>% 
  na.omit() %>% 
  mutate(litterfall = litter*prop_leaf/dt/0.45*3.65) %>% # only leaf from x/x/x to MgC/ha/year
  mutate(site = 'Paracou') %>% 
  select(date, dt, litterfall, site)
litter_tapajos <- vroom::vroom("data/fluxes/CD10_Litter_Tapajos_862/data/lba_km67_litter_archive.txt") %>% 
  mutate(date = lubridate::as_date(as.character(YYYYMMDD), format = "%Y%m%d")) %>% 
  mutate(site = 'Tapajos', litterfall = leaf, dt = dT) %>% 
  select(date, dt, litterfall, site)
bind_rows(litter_paracou, litter_tapajos) %>% 
  write_tsv(file = "outputs/litter.tsv")
Code 
obs_metrics <- read_tsv("outputs/litter.tsv") %>% 
  group_by(site, year = year(date)) %>% 
  summarise(pike = mean(c(litterfall[which.max(litterfall)],
                          litterfall[which.max(litterfall)-1],
                          litterfall[which.max(litterfall)+1]), na.rm = TRUE),
            pike_date = date[which.max(litterfall)],
            basal = mean(as.numeric(month(date) %in% 1:4)*litterfall, na.rm = TRUE)) %>%  
  mutate(pike_day = yday(pike_date), ratio = pike / basal) %>% 
  filter(basal > 0, pike_day > 100)
Code 
read_tsv("outputs/litter.tsv") %>% 
  ggplot(aes(date, litterfall, col = dt)) +
  geom_line(col = "lightgrey", alpha = 0.5) +
  geom_point() +
  theme_bw() +
  facet_wrap(~ site, scales = "free", nrow = 2) +
  scale_color_viridis_c("Time lapse") +
  xlab("") + ylab("Observed litterfall (MgC/ha/yr)") +
  geom_vline(aes(xintercept = pike_date), 
             data = obs_metrics, col = "red", linetype = "dashed") +
  geom_point(aes(x = pike_date, y = pike), 
             data = obs_metrics, col = "red", size = 2)

Fortnightly measurement of litterfall at Tapajos and Paracou. Time lapse between two sampling is over 15-days for few measurement indicated by point colour. Annual pike of litterfall are indicated by the red lines and thier value by the red points.
Code 
read_tsv("outputs/litter.tsv") %>% 
  ggplot(aes(yday(date), litterfall, col = year(date))) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 20, raw = TRUE)) +
  geom_line(aes(group = as.factor(year(date))), col = "lightgrey", alpha = 0.5) +
  geom_point() +
  theme_bw() +
  facet_wrap(~ site, scales = "free", nrow = 2) +
  scale_color_viridis_c("Year") +
  xlab("") + ylab("Observed litterfall (MgC/ha/yr)") +
  scale_x_continuous(breaks = yday(as_date(paste0("2001-", 1:12, "-1"))),
                   labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D")) +
  geom_point(aes(x = pike_day, y = pike), 
             data = obs_metrics, col = "red", size = 2)

Annual cycle of fortnightly measurement of litterfall at Tapajos and Paracou. Annual pike of litterfall are indicated by the red point.
Code 
obs_metrics %>% 
  ggplot(aes(pike_day, ratio, col = site, size = year)) +
  geom_point() +
  theme_bw() +
  ylab("Observed litterfall ratio") +
  xlab("Observed litterfall pike day")

Relations betweens Litterfall pike date and ratioi across years at Paracou and Tapajos.

Canopy structure

Canopy structure and dynamics was assessed using Leaf Area Index (LAI), Leaf Area Density (LAD) and Plant Area Density (PAD) from local observation using terrestrial or drone LIDAR remotely-sensed by MODIS. LIDAR data includes drone data from Paracou provided by Nicolas Barbier Bonal (G. Vincent & N. Barbier pers. com.) and terrestrial data from Tapajos etrieved from Smith et al. (2019): https://nph.onlinelibrary.wiley.com/doi/10.1111/nph.15726. Remotely-sensed data from MODIS were retrieved from PLUMBER2 on Research Data Australia: https://researchdata.edu.au/plumber2-forcing-evaluation-surface-models/1656048. Leaf area index was also extracted from PhenoCam of Wu et al. (2017) by X. Chen (pers. com.). Leaf area index per leaf cohort (LAI_age) was retrieved from Yang et al. (2023), from PhenoCam of Wu et al. (2017) by X. Chen (pers. com.), and reanalaysed by X. Chen et al. (in prep).

Code 
library(sf)
library(terra)
library(ncdf4)
pad <- vroom::vroom("data/fluxes/PAI_ULS_TimeSeries.txt") %>% 
  mutate(site = "Paracou", value = mean*0.68, variable = "PAI") %>% 
  select(site, date, value, variable)
dates <- readxl::read_xlsx("~/Documents/data/Tapajos/nph15726-sup-0002-notess4.xlsx", 
                  "abs_height", col_names = FALSE, skip = 0, n_max = 1) %>% 
  select(-`...1`) %>% 
  gather(key, date) %>% 
  mutate(key = as.character(1:n()))
lad <- readxl::read_xlsx("~/Documents/data/Tapajos/nph15726-sup-0002-notess4.xlsx", 
                  "abs_height", col_names = c("height", dates$key), skip = 1) %>% 
  gather(key, lad, -height) %>% 
  left_join(dates) %>% 
  mutate(date = as_date(date))
lad <- lad %>%
  group_by(date) %>% 
  summarise(value = sum(lad)) %>% 
  mutate(site = "Tapajos", variable = "LAD")
r1 <- nc_open("~/Documents/data/PLUMBER2/GF-Guy_2004-2014_FLUXNET2015_Met.nc")
r2 <- nc_open("~/Documents/data/PLUMBER2/BR-Sa3_2001-2003_FLUXNET2015_Met.nc")
lai <- list(Paracou = data.frame(lai = ncvar_get(r1, "LAI"),
             time  = as_datetime(ncatt_get(r1, "time")$units) +
             ncvar_get(r1, "time")),
     Tapajos = data.frame(lai = ncvar_get(r2, "LAI_alternative"),
             time  = as_datetime(ncatt_get(r2, "time")$units) +
             ncvar_get(r2, "time"))) %>% 
  bind_rows(.id = "site") %>% 
  group_by(site, date = date(time), variable = "LAI") %>% 
  summarise(value = mean(lai))
bind_rows(pad, lad, lai) %>% 
  write_tsv(file = "outputs/lai.tsv")
Code 
read_tsv("outputs/lai.tsv") %>% 
  bind_rows(readxl::read_xlsx("data/fluxes/Data_LAI_new.xlsx") %>% 
              mutate(date = as_date(paste0("2010-", month, "-15"))) %>% 
              mutate(site = "Tapajos", variable = "PhenoCam", 
                     value = `Camera_LAI (m^2 m^-2)`) %>% 
              select(date, value, site, variable)) %>% 
  mutate(week = week(date)) %>%  
  ggplot(aes(week, value, col = site, fill = site)) +
  geom_line(aes(group = year(date))) +
  geom_smooth(se = FALSE) +
  theme_bw() +
  theme(axis.title = element_blank(), legend.position = "bottom") +
  scale_fill_discrete("") +
  scale_color_discrete("") +
  facet_wrap(~ variable) +
  scale_x_continuous(breaks = week(as_date(paste("2000-", 1:12, "-01"))),
                     labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"))

Annual cycle of leaf area density variation at Paracou and Tapajos measured every 8-days by satellite (LAI) across ten years, every 15-days by drone in Paracou across three years, or every 15-days at Tapajos from the ground in Tapajos for three years or from PhenoCam in Tapajos for one year (Wu et al. 2016).
Code 
cols_laid <- c("young" = "#2bd22c", "mature" = "#197b14", "old" = "#ffa006")
read_tsv("data/fluxes/lai_age.tsv") %>% 
  ggplot(aes(week, m)) +
  geom_rect(aes(xmin = week, xmax = end), fill  = "#fff4e0",
            ymin = -Inf, ymax = Inf, alpha = 0.5,
            data = data.frame(week = week(as_date(c("2000-8-1", "2000-06-15"))),
                              m = -Inf,
                              end = week(as_date(c("2000-12-1", "2000-11-1"))),
                              site = c("Paracou", "Tapajos"))) +
  geom_ribbon(aes(ymin = l, ymax = h, fill = age), col = NA, alpha = 0.2) +
  geom_smooth(aes(col = age), se = FALSE) +
  theme_bw() +
  theme(legend.position = "bottom") +
  scale_fill_discrete("") +
  scale_color_discrete("") +
  facet_grid(type ~ site, scales = "free_x") +
  scale_x_continuous(breaks = week(as_date(paste("2000-", 1:12, "-01"))),
                     labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D")) + 
  ylab(expression("Leaf Area Index ["~m^2~m^{-~2}~"]")) + xlab("") +
  scale_color_manual("", values = as.vector(cols_laid[c("mature", "old", "young")])) +
  scale_fill_manual("", values = as.vector(cols_laid[c("mature", "old", "young")]))

Mean annual cycle from fortnightly means of leaf area index for Paracou and Tapajos with leaf cohorts of yound, mature and old leaves from Yang et al. (2023). Bands are the intervals of means across years, and rectangles in the background correspond to the site’s climatological dry season.
Code 
lai_age_wu <- readxl::read_xlsx("data/fluxes/Data_LAI_new.xlsx", 2) %>% 
  rename(young = `Obs_LAIyoung (m^2 m^-2)`,
         mature = `Obs_LAImature (m^2 m^-2)`,
         old = `Obs_LAIold (m^2 m^-2)`) %>% 
  mutate(site = "Tapajos", type = "Wu et al. (2016)") %>% 
  mutate(date = as_date(paste0("2010-", month, "-15"))) %>% 
  mutate(week = week(date)) %>% 
  select(site, week, type, young, mature, old) %>% 
  gather(age, m, -site, -week, -type)
lai_age_chen <- readxl::read_xlsx("data/fluxes/Data_LAI_new.xlsx", 3) %>% 
  mutate(site = recode(Site, "K67" = "Tapajos")) %>% 
  rename(young = `LAIyoung (m^2 m^-2)`,
         mature = `LAImature (m^2 m^-2)`,
         old = `LAIold (m^2 m^-2)`) %>% 
  mutate(type = "Chen et al. (in prep)") %>% 
  mutate(date = as_date(paste0("2010-", month, "-15"))) %>% 
  mutate(week = week(date)) %>% 
  select(site, week, type, young, mature, old) %>% 
  gather(age, m, -site, -week, -type)
cols_laid <- c("young" = "#2bd22c", "mature" = "#197b14", "old" = "#ffa006")
bind_rows(lai_age_wu, lai_age_chen) %>% 
  ggplot(aes(week, m)) +
  geom_rect(aes(xmin = week, xmax = end), fill  = "#fff4e0",
            ymin = -Inf, ymax = Inf, alpha = 0.5,
            data = data.frame(week = week(as_date(c("2000-8-1", "2000-06-15"))),
                              m = -Inf,
                              end = week(as_date(c("2000-12-1", "2000-11-1"))),
                              site = c("Paracou", "Tapajos"))) +
  geom_smooth(aes(col = age), se = FALSE) +
  theme_bw() +
  theme(legend.position = "bottom") +
  scale_fill_discrete("") +
  scale_color_discrete("") +
  facet_grid(type ~ site, scales = "free_x") +
  scale_x_continuous(breaks = week(as_date(paste("2000-", 1:12, "-01"))),
                     labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D")) + 
  ylab(expression("Leaf Area Index ["~m^2~m^{-~2}~"]")) + xlab("") +
  scale_color_manual("", values = as.vector(cols_laid[c("mature", "old", "young")])) +
  scale_fill_manual("", values = as.vector(cols_laid[c("mature", "old", "young")]))

Mean annual cycle from monthlu means of leaf area index for Paracou and Tapajos with leaf cohorts of yound, mature and old leaves from PhenoCam of Wu et al. (2017) by X. Chen (pers. com.), and reanalaysed by X. Chen et al. (in prep). Rectangles in the background correspond to the site’s climatological dry season.