Introduction

3-PG Forest (Physiological Processes Predicting Growth) was developed by Landsberg and Waring (1997) to bridge the gap between conventional, mensuration-based growth and yield, and process-based carbon balance models. The output variables it produces are of interest and relevance to forest managers.

3-PG calculates the radiant energy absorbed by forest canopies and converts it into biomass production. The efficiency of radiation conversion is modified by the effects of nutrition, soil drought (the model includes continuous calculation of water balance), atmospheric vapour pressure deficits and stand age.

Loading necessary packages.

library(sf)
library(dplyr)
library(nasapower)
library(ggplot2)
library(ggridges)
library(tidyr)
library(openxlsx)
library(r3PG)
library(knitr)

Create function to calculate vapour pressure deficit from average temperature and relative humidity.


calc_vpd <-function(tmp_ave,ur) {
  (0.6108*exp((17.27*tmp_ave)/(tmp_ave+237.3)))-(0.6108*exp((17.27*tmp_ave)/(tmp_ave+237.3)))*(ur/100)
}

Dowloanding and preparing Climate Data from NASA power for the municipality of Vitória da Conquista - Bahia - Brazil.

r3PG requires a table containing the information about monthly values for climatic data.


# Define the Lat Lon coordinates for an area of interest
coords_aoi <- c(-14.886761,-40.801220)

# Donwload climatic data from nasapower
climate_wider <- get_power(
  community = "ag",
  lonlat =coords_aoi, 
  pars = c("T2M_MIN_AVG", "T2M_MAX_AVG", "PRECTOTCORR_SUM", "RH2M",  "ALLSKY_SFC_SW_DWN"), 
  dates = c("1990", "2021"), 
  temporal_api = "monthly") 

# Transform the data for a longer format
climate_longer <- climate_wider%>%  select(-c(ANN))%>% pivot_longer(cols=c(-LON, -LAT,-PARAMETER, -YEAR),names_to = "month") %>%
  pivot_wider(names_from = "PARAMETER")%>%
  rename(tmp_max = T2M_MAX_AVG, tmp_min = T2M_MIN_AVG, prcp = PRECTOTCORR_SUM, ur = RH2M,
      srad = ALLSKY_SFC_SW_DWN,year = YEAR) %>%
  mutate( tmp_ave= (tmp_max+tmp_min)/2,vpd_day = calc_vpd(tmp_ave,ur),month = as.factor(month), frost_days=0)

# Update the month column for a numeric format
levels(climate_longer$month) <- as.factor(c("4", "8", "12", "2", "1", "7", "6", "3", "5","11", "10", "9"))

kable(summary(climate_longer))
LON LAT year month ur tmp_max tmp_min prcp srad tmp_ave vpd_day frost_days
Min. :-14.89 Min. :-40.8 Min. :1990 4 : 32 Min. :77.56 Min. : 9.64 Min. : 7.190 Min. : 31.64 Min. : 4.110 Min. : 8.415 Min. :0.1500 Min. :0
1st Qu.:-14.89 1st Qu.:-40.8 1st Qu.:1998 8 : 32 1st Qu.:81.61 1st Qu.:11.37 1st Qu.: 9.268 1st Qu.: 79.10 1st Qu.: 7.237 1st Qu.:10.286 1st Qu.:0.2136 1st Qu.:0
Median :-14.89 Median :-40.8 Median :2006 12 : 32 Median :83.12 Median :12.79 Median :10.755 Median :100.20 Median :12.800 Median :11.758 Median :0.2332 Median :0
Mean :-14.89 Mean :-40.8 Mean :2006 2 : 32 Mean :83.19 Mean :12.94 Mean :10.948 Mean :102.30 Mean :13.052 Mean :11.945 Mean :0.2343 Mean :0
3rd Qu.:-14.89 3rd Qu.:-40.8 3rd Qu.:2013 1 : 32 3rd Qu.:84.81 3rd Qu.:14.52 3rd Qu.:12.605 3rd Qu.:121.29 3rd Qu.:18.527 3rd Qu.:13.550 3rd Qu.:0.2550 3rd Qu.:0
Max. :-14.89 Max. :-40.8 Max. :2021 7 : 32 Max. :89.12 Max. :16.69 Max. :14.850 Max. :279.49 Max. :24.250 Max. :15.750 Max. :0.3162 Max. :0
NA NA NA (Other):192 NA NA NA NA NA NA NA NA

Defining paramenters for 3-PG simulation for AEC144 Growing simulations (Caldeira et al., 2020)

Attention to the units of measurement, see Reference manual for further details.

Download : data.input.xlsx


param <-read.xlsx("./data.input.xlsx",sheet= "parameters")
d_site <- read.xlsx('./data.input.xlsx', sheet = 'site')
d_species <- read.xlsx('./data.input.xlsx', sheet = 'species')
size_dist <- read.xlsx('./data.input.xlsx', sheet = 'sizeDist')

Executing 26 simulations of 7 years cycles for plantations between 1990 and 2015 in Vitória da Conquista - Bahia - Brazil.


# Define the name of the variable of interest from the 3PG output
i_var <- c('volume', 'lai',  'volume_mai',  'biom_stem', 'biom_root', 'biom_foliage')
i_lab <- c('Volume (m³/ha)','Leaf Area Index','Mean annual increment (m³/ha/yr)',
           'Stem biomass (Mg/ha)', 'Root biomass (Mg/ha)', 'Foliage biomass (Mg/ha')


# Create the results dataframe
df_result <- data.frame(date=as.Date(character()),species=character(),group = character(), variable =factor() ,value =numeric(), cycle=character())

# Loop through the years to create the simulations

years <- seq(1990,2021,1)
for (i in years) {
  #print(i)
  if (max(years)-i <6 ) {
      next
  }
# Subset the climate data for the present cycle
climate <-  prepare_climate(climate = climate_longer, from = paste0(i,'-01'), to = paste0(i+6,'-12')) 

# Update the planting dates for the present cycle
d_site$from <- paste0(i+1,'-01')
d_site$to <-  paste0(i+6,'-12')
d_species$planted <-  paste0(i,'-01')

# Run 3PG for the present cycle  
out_3PG <- run_3PG(
  site        = d_site, 
  species     = d_species, 
  climate     = climate, 
  thinning    = NULL,
  parameters  = param, 
  size_dist   = size_dist,
  settings    = list(light_model = 1, transp_model = 1, phys_model = 2, 
                     height_model = 1, correct_bias = 0, calculate_d13c = 0),
  check_input = TRUE, df_out = TRUE)

# Subset the outputs for the variables of interest
out_3PG_subset <- out_3PG %>%
      filter(variable %in% c(i_var,"age")) %>%
      mutate(variable = factor(variable, levels = c(i_var,"age")),
             cycle=as.character(paste0(i, " - ",i+6))) 

# Bind the present cycle's result to the results dataframe
df_result <- bind_rows(df_result,out_3PG_subset)

mai7 <-  tail(out_3PG[out_3PG$variable=="volume_mai","value"],1)%>% round(1)

# Plot the growth plots for the present cycle
out_3PG %>%
      filter(variable %in% c(i_var)) %>%
      mutate(variable = factor(variable, levels = c(i_var)),
             cycle=as.character(paste0(i, " - ",i+6)))%>% 
ggplot( aes(date, value))+
  geom_line( aes(color = species), size = 0.5)+
  facet_wrap( ~ variable, scales = 'free_y', ncol = 3, 
              labeller = labeller(variable = setNames(i_lab, i_var) )) +
  scale_color_brewer('', palette = 'Dark2') +
  theme_classic()+
  theme(legend.position="bottom")+
  xlab("Calendar date") + ylab('Value')+
  ggtitle(paste0("Cycle: ",i, " - ",i+6," MAI7 = ",mai7," m³/ha/yr"))
ggsave(filename = paste0("./plots/Cycle_",i, "-",i+6,".jpg"),device="jpg")

}

# Create a column representing the age
df_result <- df_result%>% mutate(age=rep(df_result[df_result$cycle=="1990 - 1996" & df_result$variable== "age",]$value, length(unique(df_result$cycle))*7)%>% as.vector()) 

# Transforming the results dataframe to a wider format
df_result_wider <- df_result%>% group_by(date,species,group, variable,cycle,age)%>% 
   pivot_wider(c(-species,-group,-age),names_from=variable,values_from=value )

Volume over time for 26 cycles simulations.

ggplot(  df_result_wider,aes(age, volume   ))+
  geom_line( aes(color = cycle), size = 0.1)+
  theme_classic()+
  #theme(legend.position="bottom")+
  xlab("Age (years)") + ylab('Volume (m³/ha)')

Mean Annual Increment over time for 26 cycles simulations.

ggplot( df_result_wider,aes(age, volume_mai   ))+
  geom_line( aes(color = cycle), size = 0.1)+
  theme_classic()+
  xlab("Age (years)") + ylab('MAI (m³/ha/year)')

Stem Biomass over time for 26 cycles simulations.


ggplot( df_result_wider,aes(age, biom_stem))+
  geom_line( aes(color = cycle), size = 0.1)+
  theme_classic()+
  xlab("Age (years)") + ylab('Stem Biomass (Mg/ha)')

Leaf Area Index over time for 26 cycles simulations.

ggplot( df_result_wider,aes(age, lai))+
  geom_line( aes(color = cycle), size = 0.1)+
  theme_classic()+
  xlab("Age (years)") + ylab('LAI')